Voyager I¶
In [3]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math as m
from scipy.optimize import fsolve, brentq, newton, least_squares, fmin
from matplotlib.lines import Line2D
from matplotlib import colors
from scipy.interpolate import griddata
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import ScalarFormatter
import time
import warnings
from tqdm import tqdm
import pickle
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set("jpl")
from poliastro.bodies import *
from poliastro.ephem import Ephem
from poliastro.frames import Planes
from poliastro.maneuver import Maneuver
from poliastro.twobody import Orbit
from poliastro.util import norm, time_range
from poliastro.plotting.static import StaticOrbitPlotter
from itertools import permutations, combinations
from random import shuffle
import random
import statistics
import seaborn as sns
from IPython.display import HTML
import matplotlib.animation as animation
from poliastro.core.iod import vallado
from poliastro.core.elements import rv2coe
import datetime
import os
from pathlib import Path
from IPython.display import HTML, display
from matplotlib.animation import PillowWriter
import io
from base64 import b64encode
def show(variable, units):
def get_var_name(var):
for name, value in globals().items():
if value is var:
return name
print(get_var_name(variable), "=", variable, units)
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
In [4]:
# Create an output directory handler
class OutputManager:
def __init__(self, base_dir="V1"):
# Initialise the output manager with a base directory
self.base_dir = Path(base_dir)
self.base_dir.mkdir(exist_ok=True)
print(f"Output directory set to: {self.base_dir.absolute()}")
# Save a figure with consistent settings
def save_file(self, filename, dpi=500, bbox_inches="tight", **kwargs):
# Ensure the filename has an extension
if not any(filename.endswith(ext) for ext in ['.png', '.jpg', '.pdf', '.svg']):
filename = f"{filename}.png"
# Create the full path
full_path = self.base_dir / filename
# Save the figure
plt.savefig(full_path, dpi=dpi, bbox_inches=bbox_inches, **kwargs)
print(f"Figure saved to: {full_path}")
return full_path
# Get the full path for a file.
def get_path(self, filename):
return self.base_dir / filename
# Create an instance of the output manager
output_figures = OutputManager("Voyager 1 - Figures")
output_pickle = OutputManager("Raw Results")
output_ani = OutputManager("Voyager 1 - Animations")
Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Voyager 1 - Figures Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Raw Results Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Voyager 1 - Animations
In [5]:
global_flag = False
In [6]:
mu_earth = Earth.k.to(u.km**3/u.s**2).value
mu_jupiter = Jupiter.k.to(u.km**3/u.s**2).value
mu_saturn = Saturn.k.to(u.km**3/u.s**2).value
mu_sun = Sun.k.to(u.km**3/u.s**2).value
mu_sun
Out[6]:
132712442099.00002
In [7]:
mu = mu_sun # Sun's gravitational parameter (km³/s²)
mu
Out[7]:
132712442099.00002
In [8]:
# Define parula colormap (MATLAB-style) with an option to reverse
def parula_map(reverse=False):
# RGB values
cm_data = [
[0.2081, 0.1663, 0.5292], [0.2116, 0.1898, 0.5777], [0.2123, 0.2138, 0.6270],
[0.2081, 0.2386, 0.6771], [0.1959, 0.2645, 0.7279], [0.1707, 0.2919, 0.7792],
[0.1253, 0.3242, 0.8303], [0.0591, 0.3598, 0.8683], [0.0117, 0.3875, 0.8820],
[0.0060, 0.4086, 0.8828], [0.0165, 0.4266, 0.8786], [0.0329, 0.4430, 0.8720],
[0.0498, 0.4586, 0.8641], [0.0629, 0.4737, 0.8554], [0.0723, 0.4887, 0.8467],
[0.0779, 0.5040, 0.8384], [0.0793, 0.5200, 0.8312], [0.0749, 0.5375, 0.8263],
[0.0641, 0.5570, 0.8240], [0.0488, 0.5772, 0.8228], [0.0343, 0.5966, 0.8199],
[0.0265, 0.6137, 0.8135], [0.0239, 0.6287, 0.8038], [0.0231, 0.6418, 0.7913],
[0.0228, 0.6535, 0.7768], [0.0267, 0.6642, 0.7607], [0.0384, 0.6743, 0.7436],
[0.0590, 0.6838, 0.7254], [0.0843, 0.6928, 0.7062], [0.1133, 0.7015, 0.6859],
[0.1453, 0.7098, 0.6646], [0.1801, 0.7177, 0.6424], [0.2178, 0.7250, 0.6193],
[0.2586, 0.7317, 0.5954], [0.3022, 0.7376, 0.5712], [0.3482, 0.7424, 0.5473],
[0.3953, 0.7459, 0.5244], [0.4420, 0.7481, 0.5033], [0.4871, 0.7491, 0.4840],
[0.5300, 0.7491, 0.4661], [0.5709, 0.7485, 0.4494], [0.6099, 0.7473, 0.4337],
[0.6473, 0.7456, 0.4188], [0.6834, 0.7435, 0.4044], [0.7184, 0.7411, 0.3905],
[0.7525, 0.7384, 0.3768], [0.7858, 0.7356, 0.3633], [0.8185, 0.7327, 0.3498],
[0.8507, 0.7299, 0.3360], [0.8824, 0.7274, 0.3217], [0.9139, 0.7258, 0.3063],
[0.9450, 0.7261, 0.2886], [0.9739, 0.7314, 0.2666], [0.9938, 0.7455, 0.2403],
[0.9990, 0.7653, 0.2164], [0.9955, 0.7861, 0.1967], [0.9880, 0.8066, 0.1794],
[0.9789, 0.8271, 0.1633], [0.9697, 0.8481, 0.1475], [0.9626, 0.8705, 0.1309],
[0.9589, 0.8949, 0.1132], [0.9598, 0.9218, 0.0948], [0.9661, 0.9514, 0.0755],
[0.9763, 0.9831, 0.0538]
]
if reverse:
cm_data = cm_data[::-1] # Reverse the list if specified
return colors.ListedColormap(cm_data)
1) Choose the epochs $T_1$ to $T_N$, corresponding to all the planets $P_1$ to $P_N$ in the mission. This is decided later in the optimisation algorithm.¶
In [10]:
#Voyager I real mission timing
scale = "tt"
date_departure = date_earth_jupiter = "1977-09-08 09:08"
date_flyby_jupiter = "1979-03-05 12:05"
#date_jupiter_saturn = "1979-04-24 07:33:03"
date_arrival = date_flyby_saturn = "1980-11-12 23:46"
1.1) Convert epochs to MJD¶
In [12]:
from astropy.time import Time
# e.g T_et = "1979-03-05 12:05:26"
def convert_ET(T_et):
# Defines time scale for astropy.time of ET, which approx equal to TT
T_tt = Time(T_et, scale='tt')
# Calculate TDB (only 1-2 milliseconds different from ET/TT)
T_tdb = T_tt.tdb
# Calculates the TAI (T_tai = T_tt - 32.184s)
T_tai = T_tt.tai
# Calculates the JD from epoch 1st January 4713 BC at 12:00 TT
T_jd = T_tt.jd
# Calculates MJD
T_mjd = T_jd - 2400000.5
# Calculate UTC using NASA's Sprice Time Subsystem
T_utc = T_tt.utc
return {"TDB": T_tdb, "TAI": T_tai, "JD": T_jd, "MJD": T_mjd, "UTC": T_utc}
scale = "mjd"
date_departure_mjd = convert_ET(date_departure)['MJD']
date_flyby_jupiter_mjd = convert_ET(date_flyby_jupiter)['MJD']
date_flyby_saturn_mjd = convert_ET(date_flyby_saturn)['MJD']
show(date_departure_mjd,"MJD")
show(date_flyby_jupiter_mjd,"MJD")
show(date_flyby_saturn_mjd,"MJD")
date_departure_mjd = 43394.38055555569 MJD date_flyby_jupiter_mjd = 43937.503472222015 MJD date_flyby_saturn_mjd = 44555.99027777789 MJD
1.2) Find the $\Delta t$ for each mission trajectory¶
In [14]:
def find_tof(start_date_days, end_date_days):
delta_t = (end_date_days - start_date_days) * (60*60*24)
return delta_t
delta_t_earth_jupiter = find_tof(date_departure_mjd, date_flyby_jupiter_mjd)
delta_t_jupiter_saturn = find_tof(date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
print("Earth-Jupiter:",delta_t_earth_jupiter,"s")
print("Jupiter-Saturn:",delta_t_jupiter_saturn,"s")
Earth-Jupiter: 46925819.999970496 s Jupiter-Saturn: 53437260.000027716 s
2) Using the known epochs $T_1$ to $T_N$, calculate their respective planetary state vectors $\{\vec{\mathbf{r}},\vec{\mathbf{V}}^{\,(\text{P})}\}$ using Algorithm 3.3.¶
In [16]:
def get_planet_state_vector(planet, date_time, scale):
# Convert the date_time in MJD to modern TDB time scale
date = Time(date_time, format='mjd', scale=scale).tdb
# Generate the planet's ephemeris at the given time
planet_ephem = Ephem.from_body(planet, date, plane=Planes.EARTH_ECLIPTIC)
# Extract the position and velocity vectors
r, v = planet_ephem.rv(date)
r = r.to(u.km)
v = v.to(u.km/u.s)
return r, v
scale = 'tt'
r_earth, v_earth = get_planet_state_vector(Earth, date_departure_mjd, scale)
r_jupiter, v_jupiter = get_planet_state_vector(Jupiter, date_flyby_jupiter_mjd, scale)
r_saturn, v_saturn = get_planet_state_vector(Saturn, date_flyby_saturn_mjd, scale)
In [17]:
print("Earth:")
display(r_earth, v_earth)
print("Jupiter:")
display(r_jupiter, v_jupiter)
print("Saturn:")
display(r_saturn, v_saturn)
Earth:
$[1.4652729 \times 10^{8},~-37326667,~-10360.24] \; \mathrm{km}$
$[6.7884394,~28.786762,~0.0010233513] \; \mathrm{\frac{km}{s}}$
Jupiter:
$[-4.8121511 \times 10^{8},~6.2707394 \times 10^{8},~8194396.5] \; \mathrm{km}$
$[-10.51368,~-7.3532518,~0.26572085] \; \mathrm{\frac{km}{s}}$
Saturn:
$[-1.4200499 \times 10^{9},~-54308565,~57383062] \; \mathrm{km}$
$[-0.14811167,~-9.6697442,~0.17492744] \; \mathrm{\frac{km}{s}}$
In [18]:
r_earth, v_earth = r_earth.value, v_earth.value
r_jupiter, v_jupiter = r_jupiter.value, v_jupiter.value
r_saturn, v_saturn = r_saturn.value, v_saturn.value
3) Now, to define the transfer trajectory we separate them into legs: $P_1\text{-}P_2$, $P_2\text{-}P_3$, ..., $P_{N-1}\text{-}P_N$. For each leg, we can use the solution to Lambert's Problem (Algorithm 3.4) to get its transfer velocity vectors. This fully defines each leg, where for instance $P_1\text{-}P_2$ leg, we have the following parameters $\big\{\vec{\mathbf{r}}_{1}, \vec{\mathbf{r}}_{2}, \vec{\mathbf{V}}^{\,(\text{P})}_1,\vec{\mathbf{V}}^{\,(\text{P})}_2, (\vec{\mathbf{v}})_{D}, (\vec{\mathbf{v}}_{T_2})_{A}, \Delta t\big\}$.¶
In [20]:
# from poliastro.core.iod import vallado
# # Lambert Solver
# def lambert(r1, r2, delta_t):
# v_d, v_a = vallado(mu_sun, r1, r2, delta_t, M=0, prograde=True, lowpath=True, numiter=1e6, rtol=1e-8)
# return v_d, v_a
import numpy as np
from scipy.optimize import brentq, newton, least_squares
# Known parameters [r1, r2, delta_t]
mu = 1.327e11 # Sun's gravitational parameter (km^3/s^2)
# Stumpff functions
def C(z):
if z > 0:
return (1 - np.cos(np.sqrt(z))) / z
elif z < 0:
return (np.cosh(np.sqrt(-z)) - 1) / (-z)
else:
return 0.5
def S(z):
if z > 0:
return (np.sqrt(z) - np.sin(np.sqrt(z))) / (z**1.5)
elif z < 0:
return (np.sinh(np.sqrt(-z)) - np.sqrt(-z)) / ((-z) ** 1.5)
else:
return 1 / 6
# Compute y parameter
def compute_y(z, r1, r2, A):
C_z = max(C(z), 1e-8) # Ensure no division by zero
return max(r1 + r2 + A * (z * S(z) - 1) / np.sqrt(C_z), 1e-8)
# Time-of-flight equation
def F_lambert(z, delta_t, r1, r2, A):
y_val = compute_y(z, r1, r2, A)
return ((y_val / C(z)) ** 1.5) * S(z) + A * np.sqrt(y_val) - np.sqrt(mu) * delta_t
# Lambert solver
def lambert(r1, r2, deltat):
# Step 1: Find the delta_theta between P1 and P2
r1_mag, r2_mag = np.linalg.norm(r1), np.linalg.norm(r2)
delta_theta = np.arccos(np.dot(r1, r2) / (r1_mag * r2_mag))
if np.cross(r1, r2)[2] < 0:
delta_theta = 2 * np.pi - delta_theta
# Step 2: Solve for the intermediate auxiliary variable z
A = np.sin(delta_theta) * np.sqrt(r1_mag * r2_mag / (1 - np.cos(delta_theta)))
def func(z_val):
return F_lambert(z_val, deltat, r1_mag, r2_mag, A)
def solver(func, z_min=-1e3, z_max=1e3, z_guess=0.01):
try:
z = brentq(lambda z: func(z), z_min, z_max)
except ValueError:
try:
z = newton(func, z_guess)
except RuntimeError:
z = least_squares(func, z_guess).x[0]
return z
z = solver(func)
y_val = compute_y(z, r1_mag, r2_mag, A)
# Step 3: Compute Lagrange coefficients f, g, gdot
f = 1 - y_val / r1_mag
g = A * np.sqrt(y_val / mu)
gdot = 1 - y_val / r2_mag
v1 = (r2 - f * r1) / g
v2 = (gdot * r2 - r1) / g
return v1, v2
v1_d, v2_a = lambert(r_earth, r_jupiter, delta_t_earth_jupiter)
v2_d, v3_a = lambert(r_jupiter, r_saturn, delta_t_jupiter_saturn)
print(v1_d, v2_a)
print(v2_d, v3_a)
[10.72074241 38.28749663 0.66273895] [-1.22415214e+01 3.46208225e+00 6.42429222e-03] [-20.97912678 -10.58575624 1.01828126] [-14.94072746 -13.42269638 0.82774958]
In [21]:
print(np.linalg.norm(v1_d), np.linalg.norm(v2_a))
print(np.linalg.norm(v2_d), np.linalg.norm(v3_a))
39.765637662752795 12.721670500594671 23.520605691058012 20.101723417599942
4) We can now work out $\Delta V_{injection}$, the initial injection burn from $P_1$ into transfer orbit,¶
In [23]:
deltaV_injection = v1_d - v_earth
deltaV_injection = np.linalg.norm(deltaV_injection)
deltaV_injection
Out[23]:
10.303632164583048
In [24]:
# Isp = 287.3 #s
# g0 = 9.80665 #m/s^2
# Vex = Isp*g0 / 1000 # km/s
# Vex, 'km/s'
In [25]:
# M_0_propulsion_module = 1123 #kg
# M_f_capacity_propulsion_module = 1039 #kg
# M_after_injection = M_0_propulsion_module*( 1 - np.exp(-deltaV_injection/Vex) )
# M_after_injection
In [26]:
# deltaV = deltaV_injection
# M0 = M_0_propulsion_module
# Mf = M_f_capacity_propulsion_module
# def rocket_eq(Vex):
# return Vex*np.log(M0/(M0 - Mf)) - deltaV
# Vex = fsolve(rocket_eq, 0.01)[0]
# Vex, 'km/s'
Flyby¶
5) Find $\Delta V_{p\,{(flyby)}}$, we must first find the hyperbolic excess velocities $\{\vec{\mathbf{v}}_{\infty_{in}},\vec{\mathbf{v}}_{\infty_{out}} \}$ at each flyby planet¶
6) Find the planet's inbound $a_{in}$ and outbound $a_{out}$ semi-major axis¶
7) Find its deflection (turning) angle $\delta$ from entering SOI to exiting SOI¶
Here, we can utilise the equation relating deflection (turning) angle $\delta$ to the outbound eccentricity $e_{out}$ of hyperbolic trajectory,¶
8) Hence, the radius of the closest approach to the planet $r_{p\,{(flyby)}}$ can be calculated¶
9) We now have all the required parameter to equate $\Delta V_{p\,{(flyby)}}$,¶
In [34]:
R_jupiter = 69911 #km
R_saturn = 58232 #km
R_jupiter_SOI = 48.2e6 #km
R_saturn_SOI = 54.5e6 #km
In [35]:
# Add must be more than R_planet but less that r_SOI <-- NEED TO ADD!!!!!
def find_deltaV_flyby(v_in, v_out, v_planet, mu_planet, R_planet, R_SOI):
vinf_in = v_in - v_planet
vinf_out = v_out - v_planet
a_in = - mu_planet/np.linalg.norm(vinf_in)**2
a_out = - mu_planet/np.linalg.norm(vinf_out)**2
deflection = np.arccos( np.dot(vinf_in,vinf_out)/( np.linalg.norm(vinf_in) * np.linalg.norm(vinf_out) ) )
e_out = 1/np.sin(deflection/2)
# def f(e_out):
# return (a_out/a_in * (e_out - 1))*np.sin( deflection - np.arcsin(1/e_out) ) - 1
# e_out = fsolve(f, 1.1)[0]
rp = a_out*(1-e_out)
if rp < R_planet and rp > R_SOI:
deltaV_flyby = 1e99
else:
vp_esc = np.sqrt(2*mu_planet/rp)
vp_flyby_in = np.sqrt( np.linalg.norm(vinf_in)**2 + (2 * mu_planet)/rp )
if vp_flyby_in > vp_esc:
deltaV_flyby = abs(np.sqrt( np.linalg.norm(vinf_out)**2 + (2 * mu_planet)/rp ) - np.sqrt( np.linalg.norm(vinf_in)**2 + (2 * mu_planet)/rp ) )
else:
deltaV_flyby = 1e99
return deltaV_flyby, rp
deltaV_jupiter, rp_jupiter = find_deltaV_flyby(v2_a,v2_d,v_jupiter,mu_jupiter,R_jupiter, R_jupiter_SOI)
deltaV_jupiter
Out[35]:
0.008783673764401811
10) Finally, to summate all the individual $\Delta V$ components of the initial injection and subsequent flybys we can equate,¶
In [37]:
deltaV_total = deltaV_injection + deltaV_jupiter
deltaV_total
Out[37]:
10.31241583834745
In [38]:
(date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
Out[38]:
(43394.38055555569, 43937.503472222015, 44555.99027777789)
11) From the calculated mission $\Delta V$, we can find the fuel mass $M_f$ used.¶
In [41]:
M0_Voyager = 2066 #kg
M0_propulsion_module = 1207 #kg
Mf_propulsion_module = 1039 #kg
M0_mission_module = 825 #kg
Mf_mission_module = 100 #kg
Mf_total_NASA = Mf_propulsion_module + Mf_mission_module
In [42]:
F = 68054 #N
t = 43 #s
g0 = 9.80665 # m/s2
# The Isp of the actual Voyager I mission - propulsion module
Isp_NASA = (F*t)/(Mf_propulsion_module*g0)
Isp_NASA, 's'
Out[42]:
(287.2009612891239, 's')
Specific impulse ($I_{sp}$) is a measure of a rocket engine's efficiency, defined as the impulse (thrust per unit weight flow of propellant) delivered per unit of propellant consumed.
My Real Simulation Data¶
In [45]:
def find_Isp(Mf, M0, deltaV):
Isp = (deltaV*1e3)/( np.log(M0/(M0 - Mf)) * g0)
return Isp
Isp_my_sim = find_Isp(Mf_total_NASA, M0_Voyager, deltaV_total)
Isp_my_sim, 's'
Out[45]:
(1312.1445693172518, 's')
The difference in $I_{sp}$ of NASA data and my_sim data is due to the inherit assumption in my model (Method of Patched Conics and Lambert Transfers) -- also the unpowered gravitational assist assumption.
Also, due to the fact NASA used a 2 stage system with the propulsion system ejected after injection into Earth-Jupiter transfer orbit. This left a 825kg mission module with around 100kg of fuel on it.
In [47]:
Vex = Isp_my_sim*g0 / 1000 # km/s
Vex, 'km/s'
Out[47]:
(12.867742540695026, 'km/s')
In [48]:
M0_Voyager = M0 = 2066 #kg
In [49]:
def compute_fuel_mass(M0, Vex, deltaV):
Mf = M0*( 1 - np.exp(-deltaV/Vex) )
return Mf
fuel_mass_total = compute_fuel_mass(M0, Vex, deltaV_total)
print(f"The mission deltaV of {deltaV_total:.4f}km/s caused {fuel_mass_total:.2f}kg of fuel used.")
The mission deltaV of 10.3124km/s caused 1139.00kg of fuel used.
Function: find_deltaV_mission()¶
In [52]:
def find_deltaV_mission(date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd):
#Step 1 - Find delta t
delta_t_earth_jupiter = find_tof(date_departure_mjd, date_flyby_jupiter_mjd) #s
delta_t_jupiter_saturn = find_tof(date_flyby_jupiter_mjd, date_flyby_saturn_mjd) #s
tof = delta_t_earth_jupiter + delta_t_jupiter_saturn #s
#Step 2 - Find planetary state vectors
scale = 'tt'
r_earth, v_earth = get_planet_state_vector(Earth, date_departure_mjd, scale)
r_jupiter, v_jupiter = get_planet_state_vector(Jupiter, date_flyby_jupiter_mjd, scale)
r_saturn, v_saturn = get_planet_state_vector(Saturn, date_flyby_saturn_mjd, scale)
r_earth, v_earth = r_earth.value, v_earth.value
r_jupiter, v_jupiter = r_jupiter.value, v_jupiter.value
r_saturn, v_saturn = r_saturn.value, v_saturn.value
#Step 3 - Define trajectory fully using Lambert Transfer
v1_d, v2_a = lambert(r_earth, r_jupiter, delta_t_earth_jupiter)
v2_d, v3_a = lambert(r_jupiter, r_saturn, delta_t_jupiter_saturn)
#Step 4 - Injection deltaV
deltaV_injection = v1_d - v_earth
deltaV_injection = np.linalg.norm(deltaV_injection)
#Step 5 - Find deltaV of flyby
deltaV_jupiter, rp_jupiter = find_deltaV_flyby(v2_a,v2_d,v_jupiter,mu_jupiter,R_jupiter,R_jupiter_SOI)
#if deltaV_jupiter == 1e99:
# tof = 1e99
#Step 6 - Find total deltaV of mission
deltaV_total = deltaV_injection + deltaV_jupiter
#Step 7 - Find the fuel mass used
fuel_mass_total = compute_fuel_mass(M0, Vex, deltaV_total)
data_dict = {
"deltaV_injection (km/s)": deltaV_injection,
"deltaV_jupiter (km/s)": deltaV_jupiter,
"rp_jupiter (km)": rp_jupiter,
"tof (days)": tof/(60*60*24)
}
return deltaV_total, fuel_mass_total, data_dict
deltaV_REAL, fuel_mass_REAL, data_dict_REAL = find_deltaV_mission(date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
deltaV_REAL, fuel_mass_REAL, data_dict_REAL
Out[52]:
(10.31241583834745,
1139.0,
{'deltaV_injection (km/s)': 10.303632164583048,
'deltaV_jupiter (km/s)': 0.008783673764401811,
'rp_jupiter (km)': 339964.07869084046,
'tof (days)': 1161.6097222222015})
In [53]:
dates_REAL = [date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd]
dates_REAL
Out[53]:
[43394.38055555569, 43937.503472222015, 44555.99027777789]
Optimisation Algorithms¶
In [56]:
def MJD_to_TT_calander_date(T_mjd):
T_TT = Time(T_mjd, format='mjd', scale='tt') #Time object: scale='tt' format='mjd'
return T_TT.iso
MJD_to_TT_calander_date(date_departure_mjd)
Out[56]:
'1977-09-08 09:08:00.000'
In [57]:
convert_ET( MJD_to_TT_calander_date(date_departure_mjd) )['UTC']
Out[57]:
<Time object: scale='utc' format='iso' value=1977-09-08 09:07:11.816>
In [58]:
change_in_time_days = 365
Brute force Algorithm¶
1.1) Find the $\Delta t$ for each mission trajectory¶
In [61]:
print("Earth-Jupiter:",delta_t_earth_jupiter/(60*60*24),"days")
print("Jupiter-Saturn:",delta_t_jupiter_saturn/(60*60*24),"days")
Earth-Jupiter: 543.1229166663252 days Jupiter-Saturn: 618.4868055558763 days
1.2) Define the range for each date allowed for algorithm¶
In [63]:
population = spacing = 50
print("No. of Iterations:",spacing**3)
No. of Iterations: 125000
In [64]:
time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
time_range_LIST
Out[64]:
array([-365. , -350.10204082, -335.20408163, -320.30612245,
-305.40816327, -290.51020408, -275.6122449 , -260.71428571,
-245.81632653, -230.91836735, -216.02040816, -201.12244898,
-186.2244898 , -171.32653061, -156.42857143, -141.53061224,
-126.63265306, -111.73469388, -96.83673469, -81.93877551,
-67.04081633, -52.14285714, -37.24489796, -22.34693878,
-7.44897959, 7.44897959, 22.34693878, 37.24489796,
52.14285714, 67.04081633, 81.93877551, 96.83673469,
111.73469388, 126.63265306, 141.53061224, 156.42857143,
171.32653061, 186.2244898 , 201.12244898, 216.02040816,
230.91836735, 245.81632653, 260.71428571, 275.6122449 ,
290.51020408, 305.40816327, 320.30612245, 335.20408163,
350.10204082, 365. ])
In [65]:
departure_time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
departure_time_range_LIST
Out[65]:
array([-365. , -350.10204082, -335.20408163, -320.30612245,
-305.40816327, -290.51020408, -275.6122449 , -260.71428571,
-245.81632653, -230.91836735, -216.02040816, -201.12244898,
-186.2244898 , -171.32653061, -156.42857143, -141.53061224,
-126.63265306, -111.73469388, -96.83673469, -81.93877551,
-67.04081633, -52.14285714, -37.24489796, -22.34693878,
-7.44897959, 7.44897959, 22.34693878, 37.24489796,
52.14285714, 67.04081633, 81.93877551, 96.83673469,
111.73469388, 126.63265306, 141.53061224, 156.42857143,
171.32653061, 186.2244898 , 201.12244898, 216.02040816,
230.91836735, 245.81632653, 260.71428571, 275.6122449 ,
290.51020408, 305.40816327, 320.30612245, 335.20408163,
350.10204082, 365. ])
In [66]:
flyby_jupiter_time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
flyby_jupiter_time_range_LIST
Out[66]:
array([-365. , -350.10204082, -335.20408163, -320.30612245,
-305.40816327, -290.51020408, -275.6122449 , -260.71428571,
-245.81632653, -230.91836735, -216.02040816, -201.12244898,
-186.2244898 , -171.32653061, -156.42857143, -141.53061224,
-126.63265306, -111.73469388, -96.83673469, -81.93877551,
-67.04081633, -52.14285714, -37.24489796, -22.34693878,
-7.44897959, 7.44897959, 22.34693878, 37.24489796,
52.14285714, 67.04081633, 81.93877551, 96.83673469,
111.73469388, 126.63265306, 141.53061224, 156.42857143,
171.32653061, 186.2244898 , 201.12244898, 216.02040816,
230.91836735, 245.81632653, 260.71428571, 275.6122449 ,
290.51020408, 305.40816327, 320.30612245, 335.20408163,
350.10204082, 365. ])
In [67]:
flyby_saturn_time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
flyby_saturn_time_range_LIST
Out[67]:
array([-365. , -350.10204082, -335.20408163, -320.30612245,
-305.40816327, -290.51020408, -275.6122449 , -260.71428571,
-245.81632653, -230.91836735, -216.02040816, -201.12244898,
-186.2244898 , -171.32653061, -156.42857143, -141.53061224,
-126.63265306, -111.73469388, -96.83673469, -81.93877551,
-67.04081633, -52.14285714, -37.24489796, -22.34693878,
-7.44897959, 7.44897959, 22.34693878, 37.24489796,
52.14285714, 67.04081633, 81.93877551, 96.83673469,
111.73469388, 126.63265306, 141.53061224, 156.42857143,
171.32653061, 186.2244898 , 201.12244898, 216.02040816,
230.91836735, 245.81632653, 260.71428571, 275.6122449 ,
290.51020408, 305.40816327, 320.30612245, 335.20408163,
350.10204082, 365. ])
In [68]:
deltaV_REAL, fuel_mass_REAL
Out[68]:
(10.31241583834745, 1139.0)
1.3) Run Algorithm¶
In [70]:
%%time
#warnings.filterwarnings('ignore')
if global_flag:
deltaV_list_BRUTE = []
time_range_list_index = []
dates_list = []
# Start timing
start_time = time.time()
for i in tqdm(range(0,spacing)):
date_departure_mjd_BRUTE = date_departure_mjd + departure_time_range_LIST[i]
for j in range(0,spacing):
date_flyby_jupiter_mjd_BRUTE = date_flyby_jupiter_mjd + flyby_jupiter_time_range_LIST[j]
for k in range(0,spacing):
date_flyby_saturn_mjd_BRUTE = date_flyby_saturn_mjd + flyby_saturn_time_range_LIST[k]
time_range_list_index.append([i,j,k])
try:
deltaV_mission = find_deltaV_mission(date_departure_mjd_BRUTE, date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE)
deltaV_list_BRUTE.append(deltaV_mission)
dates_list.append([date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE,
deltaV_mission])
except Exception:
deltaV_list_BRUTE.append([1e99, 1e99 ,{'temp': 1e99}])
# End timing
end_time = time.time()
execution_time_BRUTE_FORCE = end_time - start_time
print(f"\nExecution Time: {execution_time_BRUTE_FORCE:.4f} seconds")
CPU times: user 2 µs, sys: 0 ns, total: 2 µs Wall time: 4.77 µs
In [71]:
# Contour plot of the problem space --> look at local minima
# Holo View
#deltaV_list_BRUTE
#dates_list
#deltaV_list_BRUTE
In [72]:
if global_flag:
filename_brute = output_pickle.get_path("brute_results_V1.pkl")
results_brute = {
"deltaV_list_BRUTE": deltaV_list_BRUTE,
"time_range_list_index": time_range_list_index,
"dates_list": dates_list,
"execution_time_BRUTE_FORCE": execution_time_BRUTE_FORCE
}
with open(filename_brute, "wb") as f:
pickle.dump(results_brute, f)
In [73]:
filename_brute = output_pickle.get_path("brute_results_V1.pkl")
with open(filename_brute, "rb") as f:
data = pickle.load(f)
deltaV_list_BRUTE = data["deltaV_list_BRUTE"]
time_range_list_index = data["time_range_list_index"]
dates_list = data["dates_list"]
execution_time_BRUTE_FORCE = data["execution_time_BRUTE_FORCE"]
if not global_flag:
print(f"Execution Time: {execution_time_BRUTE_FORCE:.4f} seconds")
Execution Time: 412.9018 seconds
In [74]:
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN = min(deltaV_list_BRUTE, key=lambda x: x[1]) #min(deltaV_list_BRUTE)
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN
Out[74]:
(9.411191461145675,
1071.7476431066705,
{'deltaV_injection (km/s)': 9.40484890472681,
'deltaV_jupiter (km/s)': 0.006342556418864831,
'rp_jupiter (km)': 820160.096156483,
'tof (days)': 1519.1607426303601})
In [75]:
data_dict_BRUTE_MIN['tof (days)']
Out[75]:
1519.1607426303601
In [76]:
# Find the index of the minimum value
min_index = deltaV_list_BRUTE.index(min(deltaV_list_BRUTE, key=lambda x: x[0]))
min_index
Out[76]:
61798
In [77]:
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN = deltaV_list_BRUTE[min_index]
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN
Out[77]:
(9.411191461145675,
1071.7476431066705,
{'deltaV_injection (km/s)': 9.40484890472681,
'deltaV_jupiter (km/s)': 0.006342556418864831,
'rp_jupiter (km)': 820160.096156483,
'tof (days)': 1519.1607426303601})
In [78]:
change_in_time = time_range_list_index[min_index]
change_in_time
Out[78]:
[24, 35, 48]
In [79]:
date_departure_mjd_BRUTE_MIN = date_departure_mjd + departure_time_range_LIST[change_in_time[0]]
date_flyby_jupiter_mjd_BRUTE_MIN = date_flyby_jupiter_mjd + flyby_jupiter_time_range_LIST[change_in_time[1]]
date_flyby_saturn_mjd_BRUTE_MIN = date_flyby_saturn_mjd + flyby_saturn_time_range_LIST[change_in_time[2]]
dates_list_BRUTE_MIN = [date_departure_mjd_BRUTE_MIN, date_flyby_jupiter_mjd_BRUTE_MIN, date_flyby_saturn_mjd_BRUTE_MIN]
dates_list_BRUTE_MIN
Out[79]:
[43386.931575963856, 44093.93204365059, 44906.092318594216]
In [80]:
change_in_dates_list_BRUTE_MIN = [
(date_departure_mjd - date_departure_mjd_BRUTE_MIN),
(date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_BRUTE_MIN),
(date_flyby_saturn_mjd - date_flyby_saturn_mjd_BRUTE_MIN)
]
change_in_dates_list_BRUTE_MIN
Out[80]:
[7.448979591834359, -156.42857142857247, -350.102040816324]
In [81]:
print("Earth-Jupiter:",date_flyby_jupiter_mjd_BRUTE_MIN - date_departure_mjd_BRUTE_MIN,"days")
print("Jupiter-Saturn:",date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_jupiter_mjd_BRUTE_MIN,"days")
Earth-Jupiter: 707.000467686732 days Jupiter-Saturn: 812.1602749436279 days
In [82]:
date_departure_BRUTE_MIN = MJD_to_TT_calander_date(date_departure_mjd_BRUTE_MIN)
date_flyby_jupiter_BRUTE_MIN = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_BRUTE_MIN)
date_flyby_saturn_BRUTE_MIN = MJD_to_TT_calander_date(date_flyby_saturn_mjd_BRUTE_MIN)
print("Date of Departure (Min):", MJD_to_TT_calander_date(date_departure_mjd_BRUTE_MIN))
print("Date of Jupiter Flyby (Min):", MJD_to_TT_calander_date(date_flyby_jupiter_mjd_BRUTE_MIN))
print("Date of Saturn Flyby (Min):", MJD_to_TT_calander_date(date_flyby_saturn_mjd_BRUTE_MIN))
Date of Departure (Min): 1977-08-31 22:21:28.163 Date of Jupiter Flyby (Min): 1979-08-08 22:22:08.571 Date of Saturn Flyby (Min): 1981-10-29 02:12:56.327
In [83]:
# Extracting MJD values from the brute force results
departure_mjd = [dates_list[num][0] for num in range(0, len(dates_list))]
jupiter_flyby_mjd = [dates_list[num][1] for num in range(0, len(dates_list))]
saturn_flyby_mjd = [dates_list[num][2] for num in range(0, len(dates_list))]
deltaV_values = [dates_list[num][3][0] for num in range(0, len(dates_list))]
departure_mjd_change_in_time = [(element - date_departure_mjd) for element in departure_mjd]
jupiter_flyby_mjd_change_in_time = [(element - date_flyby_jupiter_mjd) for element in jupiter_flyby_mjd]
saturn_flyby_mjd_change_in_time = [(element - date_flyby_saturn_mjd) for element in saturn_flyby_mjd]
1.4) Contour Graph of Problem Space¶
In [85]:
# Filter out any extreme values
max_deltav_threshold = 20 # km/s (adjust as needed)
valid_indices = [i for i in range(len(deltaV_values)) if deltaV_values[i] < max_deltav_threshold]
# Extract valid data
valid_departure_mjd = [departure_mjd_change_in_time[i] for i in valid_indices]
valid_jupiter_mjd = [jupiter_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_saturn_mjd = [saturn_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_deltaV = [deltaV_values[i] for i in valid_indices]
# Create a 3D plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
scatter = ax.scatter(valid_departure_mjd, valid_jupiter_mjd, valid_saturn_mjd,
c=valid_deltaV, cmap=parula_map(True), s=5, alpha=0.25,
norm=colors.Normalize(vmin=min(valid_deltaV), vmax=min(valid_deltaV) * 1.5))
# Highlight the REAL delta-V point
ax.scatter([0],
[0],
[0],
color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
# Highlight the Minimum Brute Force delta-V point
ax.scatter([date_departure_mjd_BRUTE_MIN - date_departure_mjd],
[date_flyby_jupiter_mjd_BRUTE_MIN - date_flyby_jupiter_mjd],
[date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_saturn_mjd],
color='fuchsia', s=200, marker='*', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Minimum Delta-V: {deltaV_BRUTE_MIN:.2f} km/s\nMinimum Fuel Mass: {fuel_mass_BRUTE_MIN:.2f} kg')
cbar = plt.colorbar(scatter, ax=ax, label='Delta-V (km/s)')
ax.set_xlabel('Change in Earth Departure Epoch (days)')
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)')
ax.set_zlabel('Change in Saturn Flyby Epoch (days)')
ax.set_title(r"3D Visualisation of Brute Force Problem Space ($\pm$ 365 days bounds)")
limit = 375
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
ax.legend(loc = "best")
plt.tight_layout()
output_figures.save_file("voyager_I_brute_force_contour_3D.png", dpi=500, bbox_inches="tight")
#ax.view_init(elev=30, azim=45) # Adjust elevation & azimuth for a better view
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_brute_force_contour_3D.png
In [86]:
valid_departure_mjd = np.array(departure_mjd_change_in_time)
valid_jupiter_mjd = np.array(jupiter_flyby_mjd_change_in_time)
valid_saturn_mjd = np.array(saturn_flyby_mjd_change_in_time)
limit = 365
plt.close('all') # Close any existing plots
fig = plt.figure(figsize=(15, 12), constrained_layout=False)
plt.subplots_adjust(top=0.9, bottom=0.05) # Adjust top and bottom margins
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
ax = fig.add_subplot(111, projection='3d')
# X-Y Plane (Blue) Scatter and Surface
ax.scatter(valid_departure_mjd, valid_jupiter_mjd, -limit, color='blue', s=1, alpha=0.4)
grid_x, grid_y = np.mgrid[min(valid_departure_mjd):max(valid_departure_mjd):50j,
min(valid_jupiter_mjd):max(valid_jupiter_mjd):50j]
grid_z = np.full_like(grid_x, -limit)
#ax.plot_surface(grid_x, grid_y, grid_z, color='blue', alpha=0.2)
# Y-Z Plane (Red) Scatter and Surface
ax.scatter(-limit, valid_jupiter_mjd, valid_saturn_mjd, color='red', s=1, alpha=0.4)
grid_y, grid_z = np.mgrid[min(valid_jupiter_mjd):max(valid_jupiter_mjd):50j,
min(valid_saturn_mjd):max(valid_saturn_mjd):50j]
grid_x = np.full_like(grid_y, -limit)
#ax.plot_surface(grid_x, grid_y, grid_z, color='red', alpha=0.2)
# X-Z Plane (Gold) Scatter and Surface
ax.scatter(valid_departure_mjd, limit, valid_saturn_mjd, color='gold', s=1, alpha=0.4)
grid_x, grid_z = np.mgrid[min(valid_departure_mjd):max(valid_departure_mjd):50j,
min(valid_saturn_mjd):max(valid_saturn_mjd):50j]
grid_y = np.full_like(grid_x, limit)
#ax.plot_surface(grid_x, grid_y, grid_z, color='gold', alpha=0.2)
ax.set_xlabel('Change in Earth Departure Epoch (days)', labelpad=2, fontsize=14)
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)', labelpad=2, fontsize=14)
ax.set_zlabel('Change in Saturn Flyby Epoch (days)', labelpad=2, fontsize=14)
ax.set_title('3D Visualisation of Voyager I Brute Force Problem Space ($\pm$ 365 days bounds)',
fontsize=14, pad=10, y=1.02) # Slight vertical adjustment
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
output_figures.save_file("voyager_I_brute_force_problem_space.png", dpi=500, bbox_inches="tight", pad_inches=0.2)
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_brute_force_problem_space.png
In [87]:
T_real = 0
bound = 365
sampling_list = []
K=50
DeltaT = (2*bound)/K # 14.6 #days
for k in range(0,K+1):
sampling_list.append( (T_real - bound) + k*(DeltaT))
len(sampling_list)
Out[87]:
51
1.5) Fixing a Date for Problem Space 2D Contour¶
In [89]:
%%time
def optimise_trajectory(fix_parameter, spacing=50):
# Validate input parameter
if fix_parameter not in ['departure', 'jupiter', 'saturn']:
raise ValueError("Not vaild Input.")
time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
# Lists to store results
deltaV_list = []
time_range_list_index = []
dates_list = []
# Start timing
start_time = time.time()
# Determine which loops to run based on fixed parameter
if fix_parameter == 'departure':
# Fix departure date, vary Jupiter and Saturn
for j in range(spacing):
date_flyby_jupiter_mjd_BRUTE = date_flyby_jupiter_mjd + time_range_LIST[j]
for k in range(spacing):
date_flyby_saturn_mjd_BRUTE = date_flyby_saturn_mjd + time_range_LIST[k]
time_range_list_index.append([0, j, k])
try:
deltaV_mission = find_deltaV_mission(
date_departure_mjd,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE
)
deltaV_list.append(deltaV_mission)
dates_list.append([
date_departure_mjd,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE,
deltaV_mission
])
except Exception:
deltaV_list.append([1e99, 1e99, {'temp': 1e99}])
elif fix_parameter == 'jupiter':
# Fix Jupiter flyby date, vary departure and Saturn
for i in range(spacing):
date_departure_mjd_BRUTE = date_departure_mjd + time_range_LIST[i]
for k in range(spacing):
date_flyby_saturn_mjd_BRUTE = date_flyby_saturn_mjd + time_range_LIST[k]
time_range_list_index.append([i, 0, k])
try:
deltaV_mission = find_deltaV_mission(
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_BRUTE
)
deltaV_list.append(deltaV_mission)
dates_list.append([
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_BRUTE,
deltaV_mission
])
except Exception:
deltaV_list.append([1e99, 1e99, {'temp': 1e99}])
else: # fix_parameter == 'saturn'
# Fix Saturn flyby date, vary departure and Jupiter
for i in range(spacing):
date_departure_mjd_BRUTE = date_departure_mjd + time_range_LIST[i]
for j in range(spacing):
date_flyby_jupiter_mjd_BRUTE = date_flyby_jupiter_mjd + time_range_LIST[j]
time_range_list_index.append([i, j, 0])
try:
deltaV_mission = find_deltaV_mission(
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd
)
deltaV_list.append(deltaV_mission)
dates_list.append([
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd,
deltaV_mission
])
except Exception:
deltaV_list.append([1e99, 1e99, {'temp': 1e99}])
# End timing
end_time = time.time()
execution_time = end_time - start_time
# Find minimum deltaV
min_index = deltaV_list.index(min(deltaV_list, key=lambda x: x[0]))
deltaV_BRUTE_MIN_2D, fuel_mass_BRUTE_MIN_2D, data_dict_BRUTE_MIN_2D = deltaV_list[min_index]
# Get the time changes that gave minimum deltaV
change_in_time = time_range_list_index[min_index]
# Calculate the actual dates for minimum deltaV
best_dates_BRUTE_MIN_2D = [
date_departure_mjd + time_range_LIST[change_in_time[0]],
date_flyby_jupiter_mjd + time_range_LIST[change_in_time[1]],
date_flyby_saturn_mjd + time_range_LIST[change_in_time[2]]
]
# Prepare data for plotting
plot_data = {
'dates_list': dates_list,
'time_range_LIST': time_range_LIST,
'fix_parameter': fix_parameter
}
return best_dates_BRUTE_MIN_2D, deltaV_BRUTE_MIN_2D, fuel_mass_BRUTE_MIN_2D, data_dict_BRUTE_MIN_2D, execution_time, plot_data
CPU times: user 3 µs, sys: 1e+03 ns, total: 4 µs Wall time: 5.72 µs
In [90]:
def plot_optimisation_results(plot_data, max_deltav_threshold=15):
dates_list = plot_data['dates_list']
fix_parameter = plot_data['fix_parameter']
# Determine the x and y axis parameters based on fixed parameter
if fix_parameter == 'departure':
x_label = 'Change in Jupiter Flyby Epoch (days)'
y_label = 'Change in Saturn Flyby Epoch (days)'
title = '2D Contour Graph of Problem Space with Earth Departure Fixed ($\pm$365 days)'
x_mjd_ref = date_flyby_jupiter_mjd
y_mjd_ref = date_flyby_saturn_mjd
x_index = 1
y_index = 2
elif fix_parameter == 'jupiter':
x_label = 'Change in Earth Departure Epoch (days)'
y_label = 'Change in Saturn Flyby Epoch (days)'
title = '2D Contour Graph of Problem Space with Jupiter Flyby Fixed ($\pm$365 days)'
x_mjd_ref = date_departure_mjd
y_mjd_ref = date_flyby_saturn_mjd
x_index = 0
y_index = 2
else: # fix_parameter == 'saturn'
x_label = 'Change in Earth Departure Epoch (days)'
y_label = 'Change in Jupiter Flyby Epoch (days)'
title = '2D Contour Graph of Problem Space with Saturn Flyby Fixed ($\pm$365 days)'
x_mjd_ref = date_departure_mjd
y_mjd_ref = date_flyby_jupiter_mjd
x_index = 0
y_index = 1
# Extract x, y, and deltaV values
x_mjd = [dates_list[num][x_index] for num in range(len(dates_list))]
y_mjd = [dates_list[num][y_index] for num in range(len(dates_list))]
deltaV_values = [dates_list[num][3][0] for num in range(len(dates_list))]
# Convert to change in days
x_mjd_change = [(element - x_mjd_ref) for element in x_mjd]
y_mjd_change = [(element - y_mjd_ref) for element in y_mjd]
# Filter out invalid values (too high deltaV)
valid_indices = [i for i in range(len(deltaV_values)) if deltaV_values[i] < max_deltav_threshold]
# Extract valid data
valid_x_mjd = [x_mjd_change[i] for i in valid_indices]
valid_y_mjd = [y_mjd_change[i] for i in valid_indices]
valid_deltaV = [deltaV_values[i] for i in valid_indices]
# Grid sizes for interpolation
x_range = np.linspace(min(valid_x_mjd), max(valid_x_mjd), 1000)
y_range = np.linspace(min(valid_y_mjd), max(valid_y_mjd), 1000)
# Create meshgrid
X, Y = np.meshgrid(x_range, y_range)
# Interpolate deltaV values to the grid using linear interpolation
Z = griddata((valid_x_mjd, valid_y_mjd), valid_deltaV, (X, Y), method='linear')
# Create the contour plot
fig, ax = plt.subplots(figsize=(10, 8))
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
cmap = parula_map(True)
# Define contour levels within the specified range
min_deltav = np.nanmin(Z)
levels = np.linspace(min_deltav, max_deltav_threshold, 20)
# Create filled contour plot
contour = ax.contourf(X, Y, Z, levels=levels, cmap=cmap,
norm=colors.Normalize(vmin=min_deltav, vmax=max_deltav_threshold))
# Add colorbar
cbar = plt.colorbar(contour, ax=ax, label=r'$\Delta V$ (km/s)')
cbar.ax.tick_params(labelsize=12)
# Add contour lines
contour_lines = ax.contour(X, Y, Z, levels=10, colors='white', linewidths=0.5, alpha=0.7)
ax.clabel(contour_lines, inline=True, fontsize=8, fmt='%.2f')
# Find minimum deltaV point
min_idx = valid_deltaV.index(min(valid_deltaV))
min_x = valid_x_mjd[min_idx]
min_y = valid_y_mjd[min_idx]
min_dv = valid_deltaV[min_idx]
# Plot original and minimum deltaV points
ax.scatter(0, 0, color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
label=r'Real $\Delta V$: {:.2f} km/s'.format(deltaV_REAL))
#ax.scatter(min_x, min_y, color='fuchsia', s=200, marker='*', edgecolors='black', linewidth=1.2,
# label=f'Minimum $\Delta V$: {min_dv:.2f} km/s')
ax.set_xlabel(x_label, labelpad=10)
ax.set_ylabel(y_label, labelpad=10)
ax.set_title(title, pad=15)
ax.minorticks_on()
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=10)
#plt.xlim(-365,365)
#plt.ylim(-365,365)
plt.legend(loc='upper left')
plt.tight_layout()
output_figures.save_file(f"voyager_1_brute_force_{fix_parameter}_fixed_2D.png", dpi=500, bbox_inches="tight")
plt.show()
fix_parameter="departure"
BRUTE_MIN_2D_fixed_departure = optimise_trajectory("departure", 50)
BRUTE_MIN_2D_fixed_jupiter = optimise_trajectory("jupiter", 50)
BRUTE_MIN_2D_fixed_saturn = optimise_trajectory("saturn", 50)
plot_optimisation_results(BRUTE_MIN_2D_fixed_departure[-1], 15)
plot_optimisation_results(BRUTE_MIN_2D_fixed_jupiter[-1], 15)
plot_optimisation_results(BRUTE_MIN_2D_fixed_saturn[-1], 15)
/var/folders/ts/2k06pkqj7sqgrljsl3sdbxc80000gn/T/ipykernel_6847/4231714446.py:20: RuntimeWarning: overflow encountered in cosh return (np.cosh(np.sqrt(-z)) - 1) / (-z) /var/folders/ts/2k06pkqj7sqgrljsl3sdbxc80000gn/T/ipykernel_6847/4231714446.py:29: RuntimeWarning: overflow encountered in sinh return (np.sinh(np.sqrt(-z)) - np.sqrt(-z)) / ((-z) ** 1.5) /var/folders/ts/2k06pkqj7sqgrljsl3sdbxc80000gn/T/ipykernel_6847/4231714446.py:37: RuntimeWarning: invalid value encountered in scalar divide return max(r1 + r2 + A * (z * S(z) - 1) / np.sqrt(C_z), 1e-8)
Figure saved to: Voyager 1 - Figures/voyager_1_brute_force_departure_fixed_2D.png
Figure saved to: Voyager 1 - Figures/voyager_1_brute_force_jupiter_fixed_2D.png
Figure saved to: Voyager 1 - Figures/voyager_1_brute_force_saturn_fixed_2D.png
In [92]:
dates_REAL
Out[92]:
[43394.38055555569, 43937.503472222015, 44555.99027777789]
In [93]:
# Initialise an empty list to store the iterations
iterations = []
# Initial guess
initial_guess = np.array([date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd])
def objective_function(dates):
# Ensure dates stay within ±365 days of the initial guess
date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd = dates
# Ensure dates stay within ±365 days of the initial guess
lower_bound = initial_guess - 365
upper_bound = initial_guess + 365
# Check if any value is outside bounds
if np.any(dates < lower_bound) or np.any(dates > upper_bound):
deltaV = 1e99
else:
try:
deltaV, fuel_mass, data_dict = find_deltaV_mission(date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
# Store the current iteration (dates)
iterations.append(dates.copy())
except Exception:
deltaV = 1e99 # Assign a high penalty for invalid solutions
return deltaV # Minimise Delta-V
# Start timing
start_time = time.time()
# Run optimisation using `fmin` with manual bounds enforcement
optimal_dates = fmin(objective_function, initial_guess, disp=True)
# Extract optimised dates
date_departure_mjd_FMIN, date_flyby_jupiter_mjd_FMIN, date_flyby_saturn_mjd_FMIN = optimal_dates
# Compute Delta-V and fuel mass at optimised dates
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN = find_deltaV_mission(date_departure_mjd_FMIN, date_flyby_jupiter_mjd_FMIN, date_flyby_saturn_mjd_FMIN)
# Display results
print(f"Minimum Delta-V: {deltaV_FMIN:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_FMIN:.2f}")
print(f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_FMIN:.2f}")
print(f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_FMIN:.2f}")
# Convert to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_FMIN)}")
print(f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN)}")
print(f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN)}")
# End timing
end_time = time.time()
execution_time_FMIN = end_time - start_time
print(f"\nExecution Time: {execution_time_FMIN:.4f} seconds")
Optimization terminated successfully.
Current function value: 9.388030
Iterations: 201
Function evaluations: 363
Minimum Delta-V: 9.3880 km/s
Optimal Departure Date (MJD): 43388.37
Optimal Jupiter Flyby Date (MJD): 44099.45
Optimal Saturn Flyby Date (MJD): 44920.99
In calendar dates:
Departure: 1977-09-02 08:54:49.416
Jupiter flyby: 1979-08-14 10:52:57.591
Saturn flyby: 1981-11-12 23:45:59.999
Execution Time: 1.1725 seconds
In [94]:
#iterations
In [95]:
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN = find_deltaV_mission(date_departure_mjd_FMIN, date_flyby_jupiter_mjd_FMIN, date_flyby_saturn_mjd_FMIN)
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN
Out[95]:
(9.388030073904538,
1069.9564197303307,
{'deltaV_injection (km/s)': 9.388030073866002,
'deltaV_jupiter (km/s)': 3.8536285273949034e-11,
'rp_jupiter (km)': 848168.1721327668,
'tof (days)': 1532.6188724847307})
In [96]:
print("Earth-Jupiter:",date_flyby_jupiter_mjd_FMIN - date_departure_mjd_FMIN,"days")
print("Jupiter-Saturn:",date_flyby_saturn_mjd_FMIN - date_flyby_jupiter_mjd_FMIN,"days")
Earth-Jupiter: 711.0820390666777 days Jupiter-Saturn: 821.536833418053 days
In [97]:
dates_list_FMIN = [date_departure_mjd_FMIN, date_flyby_jupiter_mjd_FMIN, date_flyby_saturn_mjd_FMIN]
dates_list_FMIN
Out[97]:
[43388.371405276324, 44099.453444343, 44920.990277761055]
In [98]:
change_in_dates_list_FMIN = [
(date_departure_mjd - date_departure_mjd_FMIN),
(date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_FMIN),
(date_flyby_saturn_mjd - date_flyby_saturn_mjd_FMIN)
]
change_in_dates_list_FMIN
Out[98]:
[6.009150279365713, -161.94997212098679, -364.99999998316343]
In [99]:
date_departure_FMIN = MJD_to_TT_calander_date(date_departure_mjd_FMIN)
date_flyby_jupiter_FMIN = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN)
date_flyby_saturn_FMIN = MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN)
print("Date of Departure (Min):", MJD_to_TT_calander_date(date_departure_mjd_FMIN))
print("Date of Jupiter Flyby (Min):", MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN))
print("Date of Saturn Flyby (Min):", MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN))
Date of Departure (Min): 1977-09-02 08:54:49.416 Date of Jupiter Flyby (Min): 1979-08-14 10:52:57.591 Date of Saturn Flyby (Min): 1981-11-12 23:45:59.999
Contour Graph of Problem Space¶
In [101]:
# Extract the dates from the iterations list
departure_dates_FMIN = np.array([iter[0] for iter in iterations])
jupiter_flyby_dates_FMIN = np.array([iter[1] for iter in iterations])
saturn_flyby_dates_FMIN = np.array([iter[2] for iter in iterations])
# Calculate the corresponding deltaV values for the iterations
deltaV_values_FMIN = []
for dates in iterations:
try:
deltaV, _, _ = find_deltaV_mission(*dates)
except Exception:
deltaV_values_FMIN.append(1e99)
continue
deltaV_values_FMIN.append(deltaV)
# Create a 3D plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
# Scatter plot with color based on deltaV
scatter = ax.scatter(departure_dates_FMIN - date_departure_mjd, jupiter_flyby_dates_FMIN - date_flyby_jupiter_mjd, saturn_flyby_dates_FMIN - date_flyby_saturn_mjd,
c=deltaV_values_FMIN, cmap=parula_map(True), s=10, alpha=0.7,
norm=plt.Normalize(vmin=min(deltaV_values_FMIN), vmax=min(deltaV_values_FMIN) * 1.5))
#ax.plot(departure_dates_FMIN - date_departure_mjd, jupiter_flyby_dates_FMIN - date_flyby_jupiter_mjd, saturn_flyby_dates_FMIN - date_flyby_saturn_mjd,
# linestyle='-', color='g', markersize=6, linewidth=0.5, alpha=0.5)
ax.scatter([0],
[0],
[0],
color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
ax.scatter(date_departure_mjd_FMIN - date_departure_mjd,
date_flyby_jupiter_mjd_FMIN - date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_FMIN - date_flyby_saturn_mjd,
color='fuchsia', s=200, marker='*', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Minimum Delta-V: {deltaV_FMIN:.4f} km/s\nMinimum Fuel Mass: {fuel_mass_FMIN:.2f} kg')
cbar = plt.colorbar(scatter, ax=ax, label='Delta-V (km/s)')
ax.set_xlabel('Change in Earth Departure Epoch (days)')
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)')
ax.set_zlabel('Change in Saturn Flyby Epoch (days)')
ax.set_title('3D Visualisation of FMIN Problem Space ($\pm$ 365 days bounds)')
ax.legend(loc = "best")
# Set axis limits
# limit = 375
# ax.set_xlim(-limit, limit)
# ax.set_ylim(-limit, limit)
# ax.set_zlim(-limit, limit)
# Show the plot
plt.tight_layout()
output_figures.save_file("voyager_I_fmin_contour_3D.png", dpi=500, bbox_inches="tight")
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_fmin_contour_3D.png
In [102]:
#find_deltaV_mission(departure_dates[-1], jupiter_flyby_dates[-1], saturn_flyby_dates[-1])
change_in_dates_list_FMIN
Out[102]:
[6.009150279365713, -161.94997212098679, -364.99999998316343]
Genetic Algorithm (GA)¶
In [104]:
min_departure_date = date_departure_mjd - change_in_time_days
max_departure_date = date_departure_mjd + change_in_time_days
min_flyby_jupiter = date_flyby_jupiter_mjd - change_in_time_days
max_flyby_jupiter = date_flyby_jupiter_mjd + change_in_time_days
min_flyby_saturn = date_flyby_saturn_mjd - change_in_time_days
max_flyby_saturn = date_flyby_saturn_mjd + change_in_time_days
Step 1: Generate Initial Population¶
In [106]:
def initial_population(n_population):
# Generates initial population within mission constraints.
population = []
population.append(dates_REAL)
# Add interval to generate population
for num in range(n_population-1):
date_departure_pop = random.uniform(min_departure_date, max_departure_date)
date_flyby_jupiter_pop = random.uniform(min_flyby_jupiter, max_flyby_jupiter)
date_flyby_saturn_pop = random.uniform(min_flyby_saturn, max_flyby_saturn)
# Ensure chronological order of mission events
date_flyby_jupiter_pop = max(date_flyby_jupiter_pop, date_departure_pop)
date_flyby_saturn_pop = max(date_flyby_saturn_pop, date_flyby_jupiter_pop)
population.append([date_departure_pop, date_flyby_jupiter_pop, date_flyby_saturn_pop])
return population
#pop = initial_population(2000)
Step 2: Evaluation [Fitness Function (Objective Function)]¶
In [108]:
def fitness(individual):
# Evaluates fitness based on deltaV_total minimisation.
try:
date_departure_mjd, date_flyby_jupiter_mjd, date_flyby_saturn_mjd = individual
deltaV, fuel_mass, data_dict = find_deltaV_mission(*individual)
except Exception:
deltaV = 1e6 #np.inf # Assign a large penalty for infeasible solutions
return 1 / (deltaV + 1e-6) # Avoid division by zero
#np.array([fitness(ind) for ind in pop])
Step 3: Sort Population by Fitness¶
In [110]:
def sort_population_by_fitness(population, fitness_scores, n_population):
sorted_indices = np.argsort(fitness_scores)[::-1] # Sort by highest fitness
population = [population[i] for i in sorted_indices[:n_population]] # Ensures length population is 'n_population'
return population
In [112]:
def elitism(population, n_elite):
elite_individuals = population[:n_elite]
return elite_individuals
Step 5: Selection¶
Tournament Selection¶
In [115]:
def tournament_selection(population, fitness_scores, tournament_size=3):
# Select random individuals for the tournament
selected_indices = np.random.choice(len(population), tournament_size, replace=False)
# Find the best individual in the tournament
tournament_fitness = [fitness_scores[i] for i in selected_indices]
winner_idx = selected_indices[np.argmax(tournament_fitness)]
return population[winner_idx]
Roulette Wheel Selection¶
In [117]:
def roulette_wheel_selection(population, fitness_values):
# Calculate the normalised selection probability for each individual
selection_prob = fitness_values / sum(fitness_values)
# Calculate cumulative sum of fitness probabilities
cum_sum = selection_prob.cumsum()
# Generate a random number between 0 and 1
r = random.random()
# Find the first index where cumulative sum exceeds random number
selected_index = np.where(cum_sum >= r)[0][0]
return population[selected_index]
pop = initial_population(20)
# Evaluate fitness function
fitness_scores = np.array([fitness(ind) for ind in pop])
# Sorts population by fitness value
pop = sort_population_by_fitness(pop, fitness_scores, 20)
roulette_wheel_selection(pop, fitness_scores)
Out[117]:
[43394.38055555569, 43937.503472222015, 44555.99027777789]
Step 6: Crossover¶
Single-Point Crossover¶
In [120]:
def single_point_crossover(parent_1, parent_2):
#Performs single-point crossover.
crossover_point = random.randint(1, 2)
offspring_1 = parent_1[:crossover_point] + parent_2[crossover_point:]
offspring_2 = parent_2[:crossover_point] + parent_1[crossover_point:]
# Ensure chronological order
offspring_1[1] = max(offspring_1[1], offspring_1[0])
offspring_1[2] = max(offspring_1[2], offspring_1[1])
offspring_2[1] = max(offspring_2[1], offspring_2[0])
offspring_2[2] = max(offspring_2[2], offspring_2[1])
return offspring_1, offspring_2
Simulated Binary Crossover (SBX)¶
In [122]:
def sbx_crossover(parent_1, parent_2, eta=1.5):
offspring_1 = parent_1[:]
offspring_2 = parent_2[:]
for i in range(len(parent_1)):
u = random.random()
if u <= 0.5:
beta = (2 * u) ** (1 / (eta + 1))
else:
beta = (1 / (2 * (1 - u))) ** (1 / (eta + 1))
# Generate offspring
offspring_1[i] = 0.5 * ((1 + beta) * parent_1[i] + (1 - beta) * parent_2[i])
offspring_2[i] = 0.5 * ((1 - beta) * parent_1[i] + (1 + beta) * parent_2[i])
# Ensure chronological order (if values represent time-based sequences)
offspring_1[0] = min(offspring_1[1], offspring_1[0])
offspring_1[1] = max(offspring_1[1], offspring_1[0])
offspring_1[2] = max(offspring_1[2], offspring_1[1])
offspring_2[0] = min(offspring_2[1], offspring_2[0])
offspring_2[1] = max(offspring_2[1], offspring_2[0])
offspring_2[2] = max(offspring_2[2], offspring_2[1])
return offspring_1, offspring_2
Step 7: Mutation (Random Perturbations)¶
In [124]:
def uniform_mutation(individual):
# Create a copy to avoid modifying the original
mutated = individual.copy()
# Define mutation strength
mutation_strength = 30 # Days to perturb randomly
# Select random index to mutate
index = random.randint(0, 2)
# Apply uniform mutation
mutated[index] += random.uniform(-mutation_strength, mutation_strength)
# Ensure sequential order
mutated[0] = min(mutated[1], mutated[0])
mutated[1] = max(mutated[1], mutated[0])
mutated[2] = max(mutated[2], mutated[1])
return mutated
In [125]:
def gaussian_mutation(individual):
# Create a copy to avoid modifying the original
mutated = individual.copy()
# Define bounds for each parameter
bounds = [
(min_departure_date, max_departure_date),
(min_flyby_jupiter, max_flyby_jupiter),
(min_flyby_saturn, max_flyby_saturn)
]
# Select random index to mutate
index = random.randint(0, 2)
# Calculate parameter range
param_range = bounds[index][1] - bounds[index][0]
# Apply Gaussian mutation with standard deviation proportional to parameter range
sigma = 0.05 * param_range # 5% of the parameter range
delta = random.gauss(0, sigma)
mutated[index] += delta
# Ensure value stays within bounds
mutated[index] = max(min(mutated[index], bounds[index][1]), bounds[index][0])
# Ensure sequential order
mutated[0] = min(mutated[1], mutated[0])
mutated[1] = max(mutated[1], mutated[0])
mutated[2] = max(mutated[2], mutated[1])
return mutated
Step 8: Running GA¶
- Initialise population
- Sort population
- Start Loop
- Evaluate fitness
- Select parents (roulette wheel)
- Crossover (SBX or copy)
- Mutation
- Combine Offspring & Elite Individuals
- Filter for constraints
- Sort & truncate population
- Elitism (store best individuals)
- Does it satisfy termination condition
- If NO, then repeat loop. If YES, end loop
- Store best (optimal) individual solution
In [127]:
def run_genetic_algorithm(n_population, n_elite, crossover_per, mutation_per, n_generations, population_per_generation):
#n_generations = n_generations + 1
population = initial_population(n_population) # Generate initial population
# Evaluate fitness function
fitness_scores_temp = np.array([fitness(ind) for ind in population])
# Sorts population by fitness value
population_temp = sort_population_by_fitness(population, fitness_scores_temp, n_population)
population_per_generation.append(population_temp)
elite_individuals = []
for generation in range(0, n_generations):
start_time = time.time()
# Evaluate fitness function
fitness_scores = np.array([fitness(ind) for ind in population])
# Sorts population by fitness value
population = sort_population_by_fitness(population, fitness_scores, n_population)
# Print progress
best_individual = population[0]
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(*best_individual)
print(f"Generation {generation}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}")
# Apply elitism - preserve the best individuals
if n_elite != 0:
elite_individuals = elitism(population, n_elite)
# Select parents for potential crossover - select enough for the whole population
parents_list = []
for i in range(n_population):
parents_list.append(tournament_selection(population, fitness_scores))
# Perform crossover with crossover probability
offspring_list = []
for i in range(0, len(parents_list) - 1, 2):
# Apply crossover based on crossover probability
if random.random() < crossover_per:
offspring_1, offspring_2 = sbx_crossover(parents_list[i], parents_list[i+1])
else:
# If no crossover, children are exact copies of parents
offspring_1, offspring_2 = parents_list[i][:], parents_list[i+1][:]
offspring_list.extend([offspring_1, offspring_2])
# Perform mutation
for i in range(len(offspring_list)):
if random.random() < mutation_per:
offspring_list[i] = gaussian_mutation(offspring_list[i])
# Combine elite individuals, parents and offspring
mixed_population = offspring_list + elite_individuals
# Filter solutions that meet the mission constraints
mixed_population = [
ind for ind in mixed_population
if np.abs(ind[0] - date_departure_mjd) <= change_in_time_days
and np.abs(ind[1] - date_flyby_jupiter_mjd) <= change_in_time_days
and np.abs(ind[2] - date_flyby_saturn_mjd) <= change_in_time_days
]
# Ensure population is not empty (if all individuals were filtered out)
if not mixed_population:
mixed_population = offspring_list + elite_individuals # Restore full population if filtering removed all
# Ensure population size remains constant (selection 'n_population' best individuals)
fitness_scores = np.array([fitness(ind) for ind in mixed_population])
population = sort_population_by_fitness(mixed_population, fitness_scores, n_population)
# Store current population
population_per_generation.append(population)
end_time = time.time()
#print(end_time - start_time)
# Return the best-found solution
best_individual = population[0]
best_list = best_deltaV, best_fuel_mass, best_data = find_deltaV_mission(*best_individual)
print(f"Generation {n_generations}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}")
print("\nFinal Optimised Solution:")
print(f"Departure Date (MJD): {best_individual[0]:.2f}")
print(f"Jupiter Flyby Date (MJD): {best_individual[1]:.2f}")
print(f"Saturn Flyby Date (MJD): {best_individual[2]:.2f}")
print(f"Total deltaV: {best_deltaV:.6f} km/s")
return best_individual, best_list, population_per_generation
In [128]:
# n_population = 2000
# n_elite = 10
# crossover_per = 0.9
# mutation_per = 0.3
# n_generations = 20
n_population = 1500
n_elite = 10
crossover_per = 0.9
mutation_per = 0.35
n_generations = 50
# n_population = 300
# n_elite = 10
# crossover_per = 0.9
# mutation_per = 0.35
# n_generations = 50
In [129]:
%%time
if global_flag:
# Start timing
start_time = time.time()
# Run the Genetic Algorithm
best_individual, best_mission_data, population_per_generation = run_genetic_algorithm(n_population, n_elite, crossover_per, mutation_per, n_generations, [])
# End timing
end_time = time.time()
execution_time_GA = end_time - start_time
print(f"\nExecution Time: {execution_time_GA:.4f} seconds\n")
CPU times: user 3 µs, sys: 1 µs, total: 4 µs Wall time: 5.01 µs
In [130]:
if global_flag:
filename_GA = output_pickle.get_path("ga_results_V1.pkl")
results_GA = {
"Best Individual": best_individual,
"Best Mission Data": best_mission_data,
"Population Per Generation": population_per_generation,
"execution_time_GA": execution_time_GA
}
with open(filename_GA, "wb") as f:
pickle.dump(results_GA, f)
In [131]:
filename_GA = output_pickle.get_path("ga_results_V1.pkl")
with open(filename_GA, "rb") as f:
data = pickle.load(f)
best_individual = data["Best Individual"]
best_mission_data = data["Best Mission Data"]
population_per_generation = data["Population Per Generation"]
execution_time_GA = data["execution_time_GA"]
if not global_flag:
for gen in range(len(population_per_generation)):
best_individual = population_per_generation[gen][0]
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(*best_individual)
print(f"Generation {gen}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}")
# Return the best-found solution
best_individual = population_per_generation[-1][0]
best_list = best_deltaV, best_fuel_mass, best_data = find_deltaV_mission(*best_individual)
print(f"Generation {n_generations}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}")
print("\nFinal Optimised Solution:")
print(f"Departure Date (MJD): {best_individual[0]:.2f}")
print(f"Jupiter Flyby Date (MJD): {best_individual[1]:.2f}")
print(f"Saturn Flyby Date (MJD): {best_individual[2]:.2f}")
print(f"Total deltaV: {best_deltaV:.6f} km/s")
print(f"\nExecution Time: {execution_time_GA:.4f} seconds\n")
Generation 0: Best Delta-V = 9.673398 km/s, Fitness = 0.103376 Generation 1: Best Delta-V = 9.673398 km/s, Fitness = 0.103376 Generation 2: Best Delta-V = 9.673398 km/s, Fitness = 0.103376 Generation 3: Best Delta-V = 9.658887 km/s, Fitness = 0.103532 Generation 4: Best Delta-V = 9.541725 km/s, Fitness = 0.104803 Generation 5: Best Delta-V = 9.458555 km/s, Fitness = 0.105724 Generation 6: Best Delta-V = 9.423885 km/s, Fitness = 0.106113 Generation 7: Best Delta-V = 9.423885 km/s, Fitness = 0.106113 Generation 8: Best Delta-V = 9.414176 km/s, Fitness = 0.106223 Generation 9: Best Delta-V = 9.392526 km/s, Fitness = 0.106468 Generation 10: Best Delta-V = 9.392526 km/s, Fitness = 0.106468 Generation 11: Best Delta-V = 9.392526 km/s, Fitness = 0.106468 Generation 12: Best Delta-V = 9.389653 km/s, Fitness = 0.106500 Generation 13: Best Delta-V = 9.389653 km/s, Fitness = 0.106500 Generation 14: Best Delta-V = 9.389653 km/s, Fitness = 0.106500 Generation 15: Best Delta-V = 9.388551 km/s, Fitness = 0.106513 Generation 16: Best Delta-V = 9.388551 km/s, Fitness = 0.106513 Generation 17: Best Delta-V = 9.388551 km/s, Fitness = 0.106513 Generation 18: Best Delta-V = 9.388551 km/s, Fitness = 0.106513 Generation 19: Best Delta-V = 9.388331 km/s, Fitness = 0.106515 Generation 20: Best Delta-V = 9.388191 km/s, Fitness = 0.106517 Generation 21: Best Delta-V = 9.388191 km/s, Fitness = 0.106517 Generation 22: Best Delta-V = 9.388191 km/s, Fitness = 0.106517 Generation 23: Best Delta-V = 9.388191 km/s, Fitness = 0.106517 Generation 24: Best Delta-V = 9.388140 km/s, Fitness = 0.106517 Generation 25: Best Delta-V = 9.388060 km/s, Fitness = 0.106518 Generation 26: Best Delta-V = 9.388060 km/s, Fitness = 0.106518 Generation 27: Best Delta-V = 9.388060 km/s, Fitness = 0.106518 Generation 28: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 29: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 30: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 31: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 32: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 33: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 34: Best Delta-V = 9.388042 km/s, Fitness = 0.106518 Generation 35: Best Delta-V = 9.388039 km/s, Fitness = 0.106519 Generation 36: Best Delta-V = 9.388038 km/s, Fitness = 0.106519 Generation 37: Best Delta-V = 9.388033 km/s, Fitness = 0.106519 Generation 38: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 39: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 40: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 41: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 42: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 43: Best Delta-V = 9.388031 km/s, Fitness = 0.106519 Generation 44: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 45: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 46: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 47: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 48: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 49: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 50: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Generation 50: Best Delta-V = 9.388030 km/s, Fitness = 0.106519 Final Optimised Solution: Departure Date (MJD): 43388.37 Jupiter Flyby Date (MJD): 44099.45 Saturn Flyby Date (MJD): 44920.99 Total deltaV: 9.388030 km/s Execution Time: 363.4540 seconds
In [132]:
best_solutions_per_generation = [element[0] for element in population_per_generation]
best_solutions_per_generation, len(best_solutions_per_generation)
Out[132]:
([[43395.24936125339, 44092.91517706043, 44880.85431078622], [43395.24936125339, 44092.91517706043, 44880.85431078622], [43395.24936125339, 44092.91517706043, 44880.85431078622], [43387.26530285188, 44031.50633424819, 44750.20799249409], [43383.92821875731, 44093.27800066935, 44889.28658652461], [43386.59838086187, 44086.44605020808, 44881.147748071206], [43390.383227149854, 44095.412559876946, 44913.24275187008], [43390.383227149854, 44095.412559876946, 44913.24275187008], [43388.18242637586, 44087.54321989983, 44892.05273861021], [43388.16004720267, 44099.95270596447, 44920.99027777789], [43388.16004720267, 44099.95270596447, 44920.99027777789], [43388.16004720267, 44099.95270596447, 44920.99027777789], [43388.43621712814, 44099.65258156039, 44920.99027777789], [43388.43621712814, 44099.65258156039, 44920.99027777789], [43388.43621712814, 44099.65258156039, 44920.99027777789], [43388.5069544984, 44099.428223188144, 44920.99027777789], [43388.5069544984, 44099.428223188144, 44920.99027777789], [43388.5069544984, 44099.428223188144, 44920.99027777789], [43388.5069544984, 44099.428223188144, 44920.99027777789], [43388.33500470211, 44099.36300614009, 44920.80595897071], [43388.21305093977, 44099.403797156214, 44920.90103293918], [43388.21305093977, 44099.403797156214, 44920.90103293918], [43388.21305093977, 44099.403797156214, 44920.90103293918], [43388.21305093977, 44099.403797156214, 44920.90103293918], [43388.279109734205, 44099.438543158794, 44920.99027777789], [43388.30055327698, 44099.448358952395, 44920.99027777789], [43388.30055327698, 44099.448358952395, 44920.99027777789], [43388.30055327698, 44099.448358952395, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.383791731445, 44099.45591207383, 44920.99027777789], [43388.408632032675, 44099.45726969802, 44920.99027777789], [43388.346084432, 44099.45064558078, 44920.99027777789], [43388.35087832446, 44099.451435719064, 44920.99027777789], [43388.368657927655, 44099.45314695144, 44920.99027777789], [43388.368657927655, 44099.45314695144, 44920.99027777789], [43388.368657927655, 44099.45314695144, 44920.99027777789], [43388.368657927655, 44099.45314695144, 44920.99027777789], [43388.368657927655, 44099.45314695144, 44920.99027777789], [43388.38018698197, 44099.45411939744, 44920.98997879902], [43388.37352437623, 44099.453495186164, 44920.98993722775], [43388.37352437623, 44099.453495186164, 44920.98993722775], [43388.37352437623, 44099.453495186164, 44920.98993722775], [43388.36794380201, 44099.45314375377, 44920.990277777884], [43388.36794380201, 44099.45314375377, 44920.990277777884], [43388.37099088016, 44099.453412763745, 44920.990260873026], [43388.369157978275, 44099.45323003389, 44920.990271185954]], 51)
In [133]:
dates_REAL
Out[133]:
[43394.38055555569, 43937.503472222015, 44555.99027777789]
GA Result Analysis¶
In [135]:
deltaV_GA, fuel_mass_GA, data_dict_GA = find_deltaV_mission(*best_individual)
deltaV_GA, fuel_mass_GA, data_dict_GA
Out[135]:
(9.388030185988406,
1069.9564284063217,
{'deltaV_injection (km/s)': 9.388030110672208,
'deltaV_jupiter (km/s)': 7.531619772294107e-08,
'rp_jupiter (km)': 848166.6383611013,
'tof (days)': 1532.6211132076787})
In [136]:
date_departure_mjd_GA = best_individual[0]
date_flyby_jupiter_mjd_GA = best_individual[1]
date_flyby_saturn_mjd_GA = best_individual[2]
dates_list_GA = best_individual
dates_list_GA
Out[136]:
[43388.369157978275, 44099.45323003389, 44920.990271185954]
In [137]:
print(f"Difference for date_departure: {date_departure_mjd - date_departure_mjd_GA}")
print(f"Difference for date_flyby_jupiter: {date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_GA}")
print(f"Difference for date_flyby_saturn: {date_flyby_saturn_mjd - date_flyby_saturn_mjd_GA}\n")
change_in_dates_list_GA = [
(date_departure_mjd - date_departure_mjd_GA),
(date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_GA),
(date_flyby_saturn_mjd - date_flyby_saturn_mjd_GA)
]
change_in_dates_list_GA
Difference for date_departure: 6.011397577414755 Difference for date_flyby_jupiter: -161.94975781187532 Difference for date_flyby_saturn: -364.99999340806244
Out[137]:
[6.011397577414755, -161.94975781187532, -364.99999340806244]
In [138]:
print("Earth-Jupiter:",date_flyby_jupiter_mjd_GA - date_departure_mjd_GA,"days")
print("Jupiter-Saturn:",date_flyby_saturn_mjd_GA - date_flyby_jupiter_mjd_GA,"days")
Earth-Jupiter: 711.0840720556153 days Jupiter-Saturn: 821.5370411520635 days
In [139]:
date_departure_GA = MJD_to_TT_calander_date(date_departure_mjd_GA)
date_flyby_jupiter_GA = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA)
date_flyby_saturn_GA = MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA)
print("Date of Departure (Min):", MJD_to_TT_calander_date(date_departure_mjd_GA))
print("Date of Jupiter Flyby (Min):", MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA))
print("Date of Saturn Flyby (Min):", MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA))
print([date_departure_GA, date_flyby_jupiter_GA, date_flyby_saturn_GA])
Date of Departure (Min): 1977-09-02 08:51:35.249 Date of Jupiter Flyby (Min): 1979-08-14 10:52:39.075 Date of Saturn Flyby (Min): 1981-11-12 23:45:59.430 ['1977-09-02 08:51:35.249', '1979-08-14 10:52:39.075', '1981-11-12 23:45:59.430']
Objective Function vs No. of Generations¶
In [141]:
generations = list(range(len(best_solutions_per_generation))) # X-axis (Generation numbers)
deltaV_values_GA = [] # Y-axis (Best Delta V per generation)
fuel_mass_values = []
for solution in best_solutions_per_generation:
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(*solution) # Compute Delta V for each best solution
deltaV_values_GA.append(best_deltaV)
fuel_mass_values.append(best_fuel_mass)
# Plot the data
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
fig, ax = plt.subplots()
plt.plot(generations[:], deltaV_values_GA[:], marker='.', linestyle='--', label="Best $\Delta V$ [km/s] per Generation")
#plt.plot(generations[:], fuel_mass_values[:], marker='.', linestyle='--', label="Best Fuel Mass $(M_f)$ [kg] per Generation")
plt.xlabel("No. of Generations")
plt.ylabel("Objective Function")
plt.title("Objective Function vs No. of Generations")
plt.legend(loc="best")
# Improve grid and ticks
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=10)
plt.tight_layout()
# Show plot
plt.show()
Fitness Function vs No. of Generations¶
In [143]:
fitness_values_deltaV = 1/np.array(deltaV_values_GA)
fitness_values_fuel_mass = 1/np.array(fuel_mass_values)
# Plot the data
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
fig, ax = plt.subplots()
#plt.plot(generations, fitness_values_deltaV, marker='.', linestyle='--', label="Fitness Function $(1/\Delta V)$ per Generation")
plt.plot(generations, fitness_values_fuel_mass, marker='.', linestyle='--', label="Fitness Function $(1/M_f)$ per Generation")
plt.xlabel("No. of Generations")
plt.ylabel("Fitness Function")
plt.title("Fitness Function vs No. of Generations")
plt.legend(loc="best")
# Improve grid and ticks
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=10)
plt.tight_layout()
# Show plot
plt.show()
In [144]:
GA_list = [date_departure_mjd_GA, date_flyby_jupiter_mjd_GA, date_flyby_saturn_mjd_GA]
GA_list_calendar = [date_departure_GA, date_flyby_jupiter_GA, date_flyby_saturn_GA]
GA_list, GA_list_calendar
Out[144]:
([43388.369157978275, 44099.45323003389, 44920.990271185954], ['1977-09-02 08:51:35.249', '1979-08-14 10:52:39.075', '1981-11-12 23:45:59.430'])
Contour Graph per Generation¶
In [146]:
# Extract best solutions per generation
best_solutions = np.array(best_solutions_per_generation)
departure_dates = best_solutions[:, 0] - date_departure_mjd
jupiter_flyby_dates = best_solutions[:, 1] - date_flyby_jupiter_mjd
saturn_flyby_dates = best_solutions[:, 2] - date_flyby_saturn_mjd
# Compute deltaV for colour mapping
deltaV_values_GA = np.array([find_deltaV_mission(*sol)[0] for sol in best_solutions])
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# Scatter plot of best solutions per generation
scatter = ax.scatter(departure_dates, jupiter_flyby_dates, saturn_flyby_dates,
c=deltaV_values_GA, cmap=parula_map(True), s=50, alpha=0.5, marker='x',
norm=plt.Normalize(vmin=min(deltaV_values_GA), vmax=max(deltaV_values_GA)))
ax.scatter([0],
[0],
[0],
color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
ax.plot([0, departure_dates[0]],
[0, jupiter_flyby_dates[0]],
[0, saturn_flyby_dates[0]],
linestyle='-', color='r', markersize=6, linewidth=0.5, alpha=0.8)
ax.plot(departure_dates, jupiter_flyby_dates, saturn_flyby_dates, linestyle='-', color='r', markersize=6, linewidth=0.5, alpha=0.8)
ax.scatter(date_departure_mjd_GA - date_departure_mjd,
date_flyby_jupiter_mjd_GA - date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_GA - date_flyby_saturn_mjd,
color='fuchsia', s=150, marker='*', edgecolors='black', linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f'Minimum Delta-V: {deltaV_GA:.4f} km/s\nMinimum Fuel Mass: {fuel_mass_GA:.2f} kg')
# Colour bar
# cbar = fig.colorbar(scatter, ax=ax, shrink=0.6, pad=0.1)
# cbar.set_label('Total Delta-V (km/s)')
cbar = plt.colorbar(scatter, ax=ax, label='$\Delta$V (km/s)')
# Set labels with increased padding
ax.set_xlabel('Change in Earth Departure Epoch (days)')
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)')
ax.set_zlabel('Change in Saturn Flyby Epoch (days)')
# Title
ax.set_title('3D Visualisation of Genetic Algorithm\'s Best Solutions per Generation ($\pm$ 365 days bounds)')
# Set axis limits
#limit = 365
#ax.set_xlim(-limit, limit)
#ax.set_ylim(-limit, limit)
#ax.set_zlim(-limit, limit)
ax.legend()
plt.tight_layout()
output_figures.save_file("voyager_I_GA_per_generation_3D.png", dpi=500, bbox_inches="tight")
# Show plot
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_GA_per_generation_3D.png
In [147]:
best_solutions_per_generation = np.array(best_solutions_per_generation)
best_solutions_per_generation[:,0]
Out[147]:
array([43395.24936125, 43395.24936125, 43395.24936125, 43387.26530285,
43383.92821876, 43386.59838086, 43390.38322715, 43390.38322715,
43388.18242638, 43388.1600472 , 43388.1600472 , 43388.1600472 ,
43388.43621713, 43388.43621713, 43388.43621713, 43388.5069545 ,
43388.5069545 , 43388.5069545 , 43388.5069545 , 43388.3350047 ,
43388.21305094, 43388.21305094, 43388.21305094, 43388.21305094,
43388.27910973, 43388.30055328, 43388.30055328, 43388.30055328,
43388.38379173, 43388.38379173, 43388.38379173, 43388.38379173,
43388.38379173, 43388.38379173, 43388.38379173, 43388.40863203,
43388.34608443, 43388.35087832, 43388.36865793, 43388.36865793,
43388.36865793, 43388.36865793, 43388.36865793, 43388.38018698,
43388.37352438, 43388.37352438, 43388.37352438, 43388.3679438 ,
43388.3679438 , 43388.37099088, 43388.36915798])
In [148]:
# Extract the population for the chosen generation
selected_population = population_per_generation[0]
# Extract dates and normalise relative to nominal mission dates
departure_dates_GA = np.array([ind[0] for ind in selected_population]) - date_departure_mjd
jupiter_flyby_dates_GA = np.array([ind[1] for ind in selected_population]) - date_flyby_jupiter_mjd
saturn_flyby_dates_GA = np.array([ind[2] for ind in selected_population]) - date_flyby_saturn_mjd
# Compute deltaV for coloring
deltaV_values_GA = np.array([find_deltaV_mission(*ind)[0] for ind in selected_population])
# Create 3D scatter plot
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
scatter = ax.scatter(departure_dates_GA, jupiter_flyby_dates_GA, saturn_flyby_dates_GA,
c=deltaV_values_GA, cmap=parula_map(True), s=10, alpha=0.4,
norm=plt.Normalize(vmin=min(deltaV_values_GA), vmax=min(deltaV_values_GA) * 1.5))
ax.scatter(best_solutions_per_generation[:-6,0] - date_departure_mjd,
best_solutions_per_generation[:-6,1] - date_flyby_jupiter_mjd,
best_solutions_per_generation[:-6,2] - date_flyby_saturn_mjd,
color='#81DB24', s=20, marker='s', edgecolors='black', linewidth=1 )
ax.scatter([0],
[0],
[0],
color='red', s=200, marker='o', edgecolors='black', linewidth=1.2,
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
ax.scatter(date_departure_mjd_GA - date_departure_mjd,
date_flyby_jupiter_mjd_GA - date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_GA - date_flyby_saturn_mjd,
color='fuchsia', s=300, marker='*', edgecolors='black', linewidth=1.2,
label=f'Minimum Delta-V: {deltaV_GA:.4f} km/s\nMinimum Fuel Mass: {fuel_mass_GA:.2f} kg')
ax.plot([0, departure_dates[0]],
[0, jupiter_flyby_dates[0]],
[0, saturn_flyby_dates[0]],
linestyle='--', color='black', markersize=6, linewidth=2, alpha=1)
ax.plot(departure_dates, jupiter_flyby_dates, saturn_flyby_dates, linestyle='--', color='black', markersize=6, linewidth=2, alpha=1)
# Set labels with increased padding
ax.set_xlabel('Change in Earth Departure Epoch (days)', labelpad=10)
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)', labelpad=20)
ax.set_zlabel('Change in Saturn Flyby Epoch (days)', labelpad=20)
# Set title with increased padding
ax.set_title(f'3D Visualisation of Genetic Algorithm\'s Initial Population Problem Space ($\pm$ 365 days bounds)')
# Format the scientific notation to be more readable
for axis in [ax.xaxis, ax.yaxis, ax.zaxis]:
formatter = ScalarFormatter(useOffset=False)
formatter.set_scientific(True)
formatter.set_powerlimits((-3, 3)) # Only use scientific notation for exponents outside this range
axis.set_major_formatter(formatter)
# Position the offset text for better visibility
ax.xaxis.get_offset_text().set_position((0, -0.1))
ax.yaxis.get_offset_text().set_position((0, -0.1))
ax.zaxis.get_offset_text().set_position((0, -0.1))
# Ensure the exponent text is horizontally aligned
ax.xaxis.get_offset_text().set_ha('left')
ax.yaxis.get_offset_text().set_ha('left')
ax.zaxis.get_offset_text().set_ha('left')
# Add colorbar with adjusted position
cbar = fig.colorbar(scatter, ax=ax, shrink=0.6, pad=0.1)
cbar.set_label('Total $\Delta V$ (km/s)')
ax.legend()
plt.tight_layout()
output_figures.save_file("voyager_I_GA_inital_pop_3D.png", dpi=500, bbox_inches="tight")
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_GA_inital_pop_3D.png
In [149]:
def plot_3d_ga_contour_generation(population_per_generation, generation_idx):
# Make sure the generation index is valid
if generation_idx < 0 or generation_idx >= len(population_per_generation):
raise ValueError(f"Generation index must be between 0 and {len(population_per_generation)-1}")
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# Get the selected population
selected_population = population_per_generation[generation_idx]
# Extract dates and normalise relative to nominal mission dates
departure_dates_GA = np.array([ind[0] for ind in selected_population]) - date_departure_mjd
jupiter_flyby_dates_GA = np.array([ind[1] for ind in selected_population]) - date_flyby_jupiter_mjd
saturn_flyby_dates_GA = np.array([ind[2] for ind in selected_population]) - date_flyby_saturn_mjd
deltaV_values_GA = np.array([find_deltaV_mission(*ind)[0] for ind in selected_population])
# Create scatter plot
scatter = ax.scatter(departure_dates_GA, jupiter_flyby_dates_GA, saturn_flyby_dates_GA,
c=deltaV_values_GA, cmap=parula_map(True), s=10, alpha=0.3,
norm=plt.Normalize(vmin=min(deltaV_values_GA), vmax=min(deltaV_values_GA) * 1.5))
# Mark the best solution for the current generation
best_individual = selected_population[0] # Best individual of the generation
ax.scatter(best_individual[0] - date_departure_mjd,
best_individual[1] - date_flyby_jupiter_mjd,
best_individual[2] - date_flyby_saturn_mjd,
color='red', s=100, marker='x', label="Best Solution")
# Mark the real solution
ax.scatter([0], [0], [0],
color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
# Set labels with increased padding
ax.set_xlabel('Change in Earth Departure Epoch (days)', labelpad=10)
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)', labelpad=20)
ax.set_zlabel('Change in Saturn Flyby Epoch (days)', labelpad=20)
# Set title with increased padding
ax.set_title(f'3D Visualisation of GA Problem Space (Generation {generation_idx})', pad=20)
# Format the scientific notation to be more readable
for axis in [ax.xaxis, ax.yaxis, ax.zaxis]:
formatter = ScalarFormatter(useOffset=False)
formatter.set_scientific(True)
formatter.set_powerlimits((-3, 3)) # Only use scientific notation for exponents outside this range
axis.set_major_formatter(formatter)
# Position the offset text for better visibility
ax.xaxis.get_offset_text().set_position((0, -0.1))
ax.yaxis.get_offset_text().set_position((0, -0.1))
ax.zaxis.get_offset_text().set_position((0, -0.1))
# Ensure the exponent text is horizontally aligned
ax.xaxis.get_offset_text().set_ha('left')
ax.yaxis.get_offset_text().set_ha('left')
ax.zaxis.get_offset_text().set_ha('left')
# Add legend
ax.legend()
# Set axis limits
limit = 365
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
# Add colorbar
cbar = fig.colorbar(scatter, ax=ax, shrink=0.6, pad=0.1)
cbar.set_label('Total $\Delta V$ (km/s)')
plt.tight_layout()
#output_figures.save_file(f"voyager_I_GA_generation_{generation_idx}.png", dpi=500, bbox_inches="tight")
plt.show()
return fig
fig = plot_3d_ga_contour_generation(population_per_generation, 10)
In [150]:
deltaV_values_for_each_gen = [find_deltaV_mission(*solution)[0] for solution in best_solutions_per_generation]
deltaV_values_for_each_gen
Out[150]:
[9.673398226512557, 9.673398226512557, 9.673398226512557, 9.658887436753394, 9.541724931289206, 9.458554967466874, 9.423885054018644, 9.423885054018644, 9.414175975790753, 9.392525795434343, 9.392525795434343, 9.392525795434343, 9.389652687163839, 9.389652687163839, 9.389652687163839, 9.388551243485637, 9.388551243485637, 9.388551243485637, 9.388551243485637, 9.38833146299216, 9.388190890444521, 9.388190890444521, 9.388190890444521, 9.388190890444521, 9.388140077046343, 9.388059693351178, 9.388059693351178, 9.388059693351178, 9.388041807870483, 9.388041807870483, 9.388041807870483, 9.388041807870483, 9.388041807870483, 9.388041807870483, 9.388041807870483, 9.388038809010995, 9.38803843842005, 9.388033358155313, 9.388030674653933, 9.388030674653933, 9.388030674653933, 9.388030674653933, 9.388030674653933, 9.388030627581857, 9.388030408862573, 9.388030408862573, 9.388030408862573, 9.388030234293744, 9.388030234293744, 9.388030194026426, 9.388030185988406]
GA Results --> Fmin¶
In [152]:
%%time
# Start timing
start_time = time.time()
# Run optimisation using `fmin` with manual bounds enforcement
optimal_dates = fmin(objective_function, GA_list, disp=True)
# Extract optimised dates
date_departure_mjd_GA_FMIN_COMBINED, date_flyby_jupiter_mjd_GA_FMIN_COMBINED, date_flyby_saturn_mjd_GA_FMIN_COMBINED = optimal_dates
# Compute Delta-V and fuel mass at optimised dates
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED = find_deltaV_mission(date_departure_mjd_GA_FMIN_COMBINED, date_flyby_jupiter_mjd_GA_FMIN_COMBINED, date_flyby_saturn_mjd_GA_FMIN_COMBINED)
# Display results
print(f"Minimum Delta-V: {deltaV_GA_FMIN_COMBINED:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_GA_FMIN_COMBINED:.2f}")
print(f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_GA_FMIN_COMBINED:.2f}")
print(f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_GA_FMIN_COMBINED:.2f}")
# Convert to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED)}")
print(f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED)}")
print(f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED)}")
# End timing
end_time = time.time()
execution_time_GA_FMIN_COMBINED = end_time - start_time
print(f"\nExecution Time: {execution_time_GA_FMIN_COMBINED:.4f} seconds")
Optimization terminated successfully.
Current function value: 9.388030
Iterations: 101
Function evaluations: 193
Minimum Delta-V: 9.3880 km/s
Optimal Departure Date (MJD): 43388.37
Optimal Jupiter Flyby Date (MJD): 44099.45
Optimal Saturn Flyby Date (MJD): 44920.99
In calendar dates:
Departure: 1977-09-02 08:51:39.410
Jupiter flyby: 1979-08-14 10:52:38.296
Saturn flyby: 1981-11-12 23:45:54.227
Execution Time: 0.5063 seconds
CPU times: user 494 ms, sys: 21.2 ms, total: 515 ms
Wall time: 506 ms
In [153]:
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED = find_deltaV_mission(date_departure_mjd_GA_FMIN_COMBINED, date_flyby_jupiter_mjd_GA_FMIN_COMBINED, date_flyby_saturn_mjd_GA_FMIN_COMBINED)
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED
Out[153]:
(9.388030168438553,
1069.9564270478534,
{'deltaV_injection (km/s)': 9.388030134667423,
'deltaV_jupiter (km/s)': 3.377112989255693e-08,
'rp_jupiter (km)': 848166.5346901253,
'tof (days)': 1532.621004831104})
In [154]:
change_in_dates_list_FMIN_COMBINED = [
(date_departure_mjd - date_departure_mjd_GA_FMIN_COMBINED),
(date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_GA_FMIN_COMBINED),
(date_flyby_saturn_mjd - date_flyby_saturn_mjd_GA_FMIN_COMBINED)
]
change_in_dates_list_FMIN_COMBINED
Out[154]:
[6.011349423162756, -161.94974879429355, -364.9999331857398]
In [155]:
print("Earth-Jupiter:",date_flyby_jupiter_mjd_GA_FMIN_COMBINED - date_departure_mjd_GA_FMIN_COMBINED,"days")
print("Jupiter-Saturn:",date_flyby_saturn_mjd_GA_FMIN_COMBINED - date_flyby_jupiter_mjd_GA_FMIN_COMBINED,"days")
Earth-Jupiter: 711.0840148837815 days Jupiter-Saturn: 821.5369899473226 days
In [156]:
date_departure_GA_FMIN_COMBINED = MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED)
date_flyby_jupiter_GA_FMIN_COMBINED = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED)
date_flyby_saturn_GA_FMIN_COMBINED = MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED)
print("Date of Departure (Min):", MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED))
print("Date of Jupiter Flyby (Min):", MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED))
print("Date of Saturn Flyby (Min):", MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED))
Date of Departure (Min): 1977-09-02 08:51:39.410 Date of Jupiter Flyby (Min): 1979-08-14 10:52:38.296 Date of Saturn Flyby (Min): 1981-11-12 23:45:54.227
Graphs¶
Comparsion Graphs¶
In [159]:
def plot_VoyagerI(date_departure, date_flyby_jupiter, date_flyby_saturn, colour, label_suffix=""):
# Convert input dates to Astropy Time in TDB scale
date_departure = Time(date_departure, scale="tt").tdb
date_flyby_jupiter = Time(date_flyby_jupiter, scale="tt").tdb
date_flyby_saturn = Time(date_flyby_saturn, scale="tt").tdb
date_arrival = date_flyby_saturn
# Define planetary ephemerides
earth = Ephem.from_body(Earth, time_range(date_departure, end=date_arrival, periods=500))
ss_earth = Orbit.from_ephem(Sun, earth, date_departure)
jupiter = Ephem.from_body(Jupiter, time_range(date_departure, end=date_arrival, periods=500))
ss_jupiter = Orbit.from_ephem(Sun, jupiter, date_flyby_jupiter)
saturn = Ephem.from_body(Saturn, time_range(date_departure, end=date_arrival, periods=500))
ss_saturn = Orbit.from_ephem(Sun, saturn, date_arrival)
# Solving for maneuver to Jupiter
man_flyby_jupiter = Maneuver.lambert(ss_earth, ss_jupiter)
ic1 = ss_earth.apply_maneuver(man_flyby_jupiter)
ic1_end = ic1.propagate(date_flyby_jupiter)
# Solving for maneuver to Saturn
man_flyby_saturn = Maneuver.lambert(ic1_end, ss_saturn)
ic2 = ic1_end.apply_maneuver(man_flyby_saturn)
ic2_end = ic2.propagate(date_arrival)
# Plot Earth's position at departure
plotter.plot_body_orbit(
Earth, date_departure, label=f"Earth Departure {label_suffix}", trail=True
)
# Plot cruise to Jupiter
plotter.plot_maneuver(
ss_earth, man_flyby_jupiter, label=f"Cruise to Jupiter {label_suffix}", color=colour
)
# Plot Jupiter's position at flyby
plotter.plot_body_orbit(
Jupiter, date_flyby_jupiter, label=f"Jupiter Flyby {label_suffix}", trail=True
)
# Plot cruise to Saturn
plotter.plot_maneuver(
ic1_end, man_flyby_saturn, label=f"Cruise to Saturn {label_suffix}", color=colour
)
# Plot Earth's position at departure
plotter.plot_body_orbit(
Earth, date_departure, label=f"Earth Departure {label_suffix}", trail=True
)
# Plot Jupiter's position at flyby
plotter.plot_body_orbit(
Jupiter, date_flyby_jupiter, label=f"Jupiter Flyby {label_suffix}", trail=True
)
# Plot Saturn's position at arrival
plotter.plot_body_orbit(
Saturn, date_flyby_saturn, label=f"Saturn Arrival {label_suffix}", trail=True
)
2-Axis¶
In [161]:
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
# Create figure and axis
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111)
plotter = StaticOrbitPlotter(ax)
# Define different line styles for each plot
colours = {
"Real Mission": "red",
"Brute Force Algorithm": "cyan",
"FMIN": "magenta",
"GA": "gold",
"GA --> FMIN": "lime"
}
# Plot missions
plot_VoyagerI(date_departure, date_flyby_jupiter, date_flyby_saturn, colours["Real Mission"], label_suffix="(Real Mission)")
plot_VoyagerI(date_departure_BRUTE_MIN, date_flyby_jupiter_BRUTE_MIN, date_flyby_saturn_BRUTE_MIN, colours["Brute Force Algorithm"], label_suffix="(Brute Force Algorithm)")
plot_VoyagerI(date_departure_FMIN, date_flyby_jupiter_FMIN, date_flyby_saturn_FMIN, colours["FMIN"], label_suffix="(FMIN)")
plot_VoyagerI(date_departure_GA, date_flyby_jupiter_GA, date_flyby_saturn_GA, colours["GA"], label_suffix="(GA)")
#plot_VoyagerI(date_departure_GA_FMIN_COMBINED, date_flyby_jupiter_GA_FMIN_COMBINED, date_flyby_saturn_GA_FMIN_COMBINED, colours["GA --> FMIN"], label_suffix="(GA --> FMIN)")
# Set fixed limits for primary axis
limit = 1.7e9
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect('equal') # Square aspect ratio
# Conversion factor (1 AU = 149,597,870.691 km)
km_to_AU = 1 / 149597870.691
# Create second set of axes for AU scale
# These won't be linked to the primary axes in the problematic way
ax_top = fig.add_axes(ax.get_position(), frameon=False)
ax_top.xaxis.tick_top()
ax_top.yaxis.tick_right()
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
ax_top.set_xlabel("x (AU)")#, labelpad=10)
ax_top.set_ylabel("y (AU)")#, labelpad=10)
ax_top.xaxis.set_label_position('top')
ax_top.yaxis.set_label_position('right')
# # Add grid
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
#ax_top.grid(which='major', color='#DDDDDD', linewidth=0.8)
#ax_top.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax_top.minorticks_on()
# Customise tick parameters
ax.tick_params(axis='both', which='major', labelsize=10)
ax.tick_params(axis='both', which='minor', labelsize=8)
ax_top.tick_params(axis='both', which='major', labelsize=10)
ax_top.tick_params(axis='both', which='minor', labelsize=8)
# Add title
plt.title(r"\textbf{Voyager I Mission Trajectory Optimisation}", pad=10)
# Function to update AU axes when primary axes change
def update_au_axes(event=None):
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
fig.canvas.draw_idle()
# Connect the update function to the 'xlim_changed' and 'ylim_changed' events
ax.callbacks.connect('xlim_changed', update_au_axes)
ax.callbacks.connect('ylim_changed', update_au_axes)
# Define different line styles for each plot
line_styles = {
"Real Mission": "-",
"Brute Force Algorithm": "-",
"FMIN": "-.",
"GA": "--",
"GA --> FMIN": ":"
}
# Define different line widths for clarity
line_widths = {
"Real Mission": 1.75,
"Brute Force Algorithm": 4,
"FMIN": 3,
"GA": 2,
"GA --> FMIN": 2
}
# Apply line styles
for line in ax.get_lines():
label = line.get_label()
for key in line_styles:
if key in label:
line.set_linestyle(line_styles[key])
line.set_linewidth(line_widths[key])
# Create custom legend
handles = [
Line2D([0], [0], color=colours['Real Mission'], ls=line_styles['Real Mission'], lw=line_widths['Real Mission'], label="Real Mission"),
Line2D([0], [0], color=colours['Brute Force Algorithm'], ls=line_styles['Brute Force Algorithm'], lw=line_widths['Brute Force Algorithm'], label="Brute Force Algorithm"),
Line2D([0], [0], color=colours['FMIN'], ls=line_styles['FMIN'], lw=line_widths['FMIN'], label="FMIN"),
Line2D([0], [0], color=colours['GA'], ls=line_styles['GA'], lw=line_widths['GA'], label="GA"),
#Line2D([0], [0], color=colours['GA --> FMIN'], ls=line_styles['GA --> FMIN'], lw=line_widths['GA --> FMIN'], label="GA \& FMIN")
]
ax.plot(0, 0, 'o', markersize=18, markerfacecolor='#FFFF00', markeredgecolor='#FF8C00', markeredgewidth=1, zorder=100)
# Add legend
ax.legend(handles=handles, loc="upper left", ncol=1, frameon=True, edgecolor='black')
# Planet colors
planet_colors = {
"Earth": "#3366FF", # Blue
"Jupiter": "#FF9933", # Orange
"Saturn": "#FFCC33" # Yellow
}
# Use the new label positions provided
planet_label_positions = {
"Earth": (-0.35e9, -0.2e9),
"Jupiter": (0.8e9, 0.5e9),
"Saturn": (0e9, 1.55e9)
}
# Add each planet marker and label
for planet, position in planet_label_positions.items():
x, y = position
# Use the absolute positions for labels instead of offsets
label_x, label_y = planet_label_positions[planet]
bbox_props = dict(boxstyle="round,pad=0.3", fc="white", ec="black", alpha=0.8)
ax.annotate(planet,
xy=(x, y), # Position of the planet
xytext=(label_x, label_y), # Absolute position for the label
textcoords='data', # Changed from 'offset points' to 'data'
bbox=bbox_props,
fontsize=11,
fontweight='bold',
zorder=101, # Ensure labels appear on top of everything
) # Add arrow connecting label to planet
# Save figure with padding
output_figures.save_file("voyager_I_trajectory_km_AU.png", dpi=500, bbox_inches="tight")
# Show plot
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_trajectory_km_AU.png
In [162]:
dates_REAL
Out[162]:
[43394.38055555569, 43937.503472222015, 44555.99027777789]
2-Axis Zoomed¶
In [164]:
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
# Create figure and axis
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111)
plotter = StaticOrbitPlotter(ax)
# Plot missions
plot_VoyagerI(date_departure, date_flyby_jupiter, date_flyby_saturn, colours["Real Mission"], label_suffix="(Real Mission)")
plot_VoyagerI(date_departure_BRUTE_MIN, date_flyby_jupiter_BRUTE_MIN, date_flyby_saturn_BRUTE_MIN, colours["Brute Force Algorithm"], label_suffix="(Brute Force Algorithm)")
plot_VoyagerI(date_departure_FMIN, date_flyby_jupiter_FMIN, date_flyby_saturn_FMIN, colours["FMIN"], label_suffix="(FMIN)")
plot_VoyagerI(date_departure_GA, date_flyby_jupiter_GA, date_flyby_saturn_GA, colours["GA"], label_suffix="(GA)")
#plot_VoyagerI(date_departure_GA_FMIN_COMBINED, date_flyby_jupiter_GA_FMIN_COMBINED, date_flyby_saturn_GA_FMIN_COMBINED, colours["GA --> FMIN"], label_suffix="(GA --> FMIN)")
# Set fixed limits for primary axis
ax.set_xlim(-.5e9, 1e9)
ax.set_ylim(-.5e9, 1.75e9)
ax.set_aspect('equal') # Square aspect ratio
# Conversion factor (1 AU = 149,597,870.691 km)
km_to_AU = 1 / 149597870.691
# Create second set of axes for AU scale
ax_top = fig.add_axes(ax.get_position(), frameon=False)
ax_top.xaxis.tick_top()
ax_top.yaxis.tick_right()
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
ax_top.set_xlabel("x (AU)")#, labelpad=10)
ax_top.set_ylabel("y (AU)")#, labelpad=10)
ax_top.xaxis.set_label_position('top')
ax_top.yaxis.set_label_position('right')
# Add grid
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
#ax_top.grid(which='major', color='#DDDDDD', linewidth=0.8)
#ax_top.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax_top.minorticks_on()
# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=10)
ax.tick_params(axis='both', which='minor', labelsize=8)
ax_top.tick_params(axis='both', which='major', labelsize=10)
ax_top.tick_params(axis='both', which='minor', labelsize=8)
# Add title
plt.title(r"\textbf{Voyager I Mission Trajectory Optimisation}", pad=10)
# Function to update AU axes when primary axes change
def update_au_axes(event=None):
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
fig.canvas.draw_idle()
# Connect the update function to the 'xlim_changed' and 'ylim_changed' events
ax.callbacks.connect('xlim_changed', update_au_axes)
ax.callbacks.connect('ylim_changed', update_au_axes)
# Apply line styles
for line in ax.get_lines():
label = line.get_label()
for key in line_styles:
if key in label:
line.set_linestyle(line_styles[key])
line.set_linewidth(line_widths[key])
# Create custom legend
handles = [
Line2D([0], [0], color=colours['Real Mission'], ls=line_styles['Real Mission'], lw=line_widths['Real Mission'], label="Real Mission"),
Line2D([0], [0], color=colours['Brute Force Algorithm'], ls=line_styles['Brute Force Algorithm'], lw=line_widths['Brute Force Algorithm'], label="Brute Force Algorithm"),
Line2D([0], [0], color=colours['FMIN'], ls=line_styles['FMIN'], lw=line_widths['FMIN'], label="FMIN"),
Line2D([0], [0], color=colours['GA'], ls=line_styles['GA'], lw=line_widths['GA'], label="GA"),
#Line2D([0], [0], color=colours['GA --> FMIN'], ls=line_styles['GA --> FMIN'], lw=line_widths['GA --> FMIN'], label="GA \& FMIN")
]
ax.plot(0, 0, 'o', markersize=18, markerfacecolor='#FFFF00', markeredgecolor='#FF8C00', markeredgewidth=1, zorder=100)
# Add legend
ax.legend(handles=handles, loc="best", ncol=1, frameon=True, edgecolor='black')
#plt.tight_layout()
output_figures.save_file("voyager_I_trajectory_zoom.png", dpi=500, bbox_inches="tight")
# Show plot
plt.show()
Figure saved to: Voyager 1 - Figures/voyager_I_trajectory_zoom.png
Genetic Algorithm (GA) Animation¶
Orbital Trajectory per Generation¶
In [167]:
# Ensure real mission is first in the best solutions list
best_solutions_per_generation,len(best_solutions_per_generation)
Out[167]:
(array([[43395.24936125, 44092.91517706, 44880.85431079],
[43395.24936125, 44092.91517706, 44880.85431079],
[43395.24936125, 44092.91517706, 44880.85431079],
[43387.26530285, 44031.50633425, 44750.20799249],
[43383.92821876, 44093.27800067, 44889.28658652],
[43386.59838086, 44086.44605021, 44881.14774807],
[43390.38322715, 44095.41255988, 44913.24275187],
[43390.38322715, 44095.41255988, 44913.24275187],
[43388.18242638, 44087.5432199 , 44892.05273861],
[43388.1600472 , 44099.95270596, 44920.99027778],
[43388.1600472 , 44099.95270596, 44920.99027778],
[43388.1600472 , 44099.95270596, 44920.99027778],
[43388.43621713, 44099.65258156, 44920.99027778],
[43388.43621713, 44099.65258156, 44920.99027778],
[43388.43621713, 44099.65258156, 44920.99027778],
[43388.5069545 , 44099.42822319, 44920.99027778],
[43388.5069545 , 44099.42822319, 44920.99027778],
[43388.5069545 , 44099.42822319, 44920.99027778],
[43388.5069545 , 44099.42822319, 44920.99027778],
[43388.3350047 , 44099.36300614, 44920.80595897],
[43388.21305094, 44099.40379716, 44920.90103294],
[43388.21305094, 44099.40379716, 44920.90103294],
[43388.21305094, 44099.40379716, 44920.90103294],
[43388.21305094, 44099.40379716, 44920.90103294],
[43388.27910973, 44099.43854316, 44920.99027778],
[43388.30055328, 44099.44835895, 44920.99027778],
[43388.30055328, 44099.44835895, 44920.99027778],
[43388.30055328, 44099.44835895, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.38379173, 44099.45591207, 44920.99027778],
[43388.40863203, 44099.4572697 , 44920.99027778],
[43388.34608443, 44099.45064558, 44920.99027778],
[43388.35087832, 44099.45143572, 44920.99027778],
[43388.36865793, 44099.45314695, 44920.99027778],
[43388.36865793, 44099.45314695, 44920.99027778],
[43388.36865793, 44099.45314695, 44920.99027778],
[43388.36865793, 44099.45314695, 44920.99027778],
[43388.36865793, 44099.45314695, 44920.99027778],
[43388.38018698, 44099.4541194 , 44920.9899788 ],
[43388.37352438, 44099.45349519, 44920.98993723],
[43388.37352438, 44099.45349519, 44920.98993723],
[43388.37352438, 44099.45349519, 44920.98993723],
[43388.3679438 , 44099.45314375, 44920.99027778],
[43388.3679438 , 44099.45314375, 44920.99027778],
[43388.37099088, 44099.45341276, 44920.99026087],
[43388.36915798, 44099.45323003, 44920.99027119]]),
51)
In [168]:
fig, ax = plt.subplots(figsize=(9, 7))
plt.close(fig) # This prevents the initial empty figure from showing
plotter = StaticOrbitPlotter(ax)
# Store previously plotted trajectories
plotted_trajectories = []
# Select a colour map (viridis, plasma, coolwarm, rainbow, etc.)
cmap = plt.cm.rainbow # Change this for different color gradients
# Flag to track if we've plotted the real mission trajectory
real_mission_plotted = False
# Function to update the plot at each generation
def update(frame):
global real_mission_plotted
# Plot the real mission trajectory only once at the first frame
if frame == 0:
# Plot real mission in red
real_traj = plot_VoyagerI(date_departure, date_flyby_jupiter, date_flyby_saturn,
"red", label_suffix="(Real Mission)")
plotted_trajectories.append(real_traj)
real_mission_plotted = True
else:
# Get best trajectory for this generation
best_traj = best_solutions_per_generation[frame-1]
date_departure_gen, date_flyby_jupiter_gen, date_flyby_saturn_gen = best_traj
# Convert MJD to calendar date
date_departure_gen = MJD_to_TT_calander_date(date_departure_gen)
date_flyby_jupiter_gen = MJD_to_TT_calander_date(date_flyby_jupiter_gen)
date_flyby_saturn_gen = MJD_to_TT_calander_date(date_flyby_saturn_gen)
# Normalise color based on generation number
norm_color = cmap(frame / len(best_solutions_per_generation))
# Plot new trajectory with gradient color
traj_plot = plot_VoyagerI(date_departure_gen, date_flyby_jupiter_gen,
date_flyby_saturn_gen, norm_color,
label_suffix=f"(Gen {frame-1})")
# Store the trajectory plot for reference
plotted_trajectories.append(traj_plot)
# Hide the legend during animation to prevent clutter
ax.legend().set_visible(False)
if frame == 0:
ax.set_title(f"Voyager I GA Optimisation - Real Mission")
else:
# Set title dynamically
ax.set_title(f"Voyager I GA Optimisation - Generation {frame-1}")
return plotted_trajectories
# Improve grid and ticks
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=10)
plt.tight_layout()
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(best_solutions_per_generation)+1,
repeat=False, blit=False)
# Save the animation as a GIF or MP4
#ani.save("voyager_I_GA_optimisation.mp4", writer="ffmpeg", fps=5)#, dpi=200, bitrate=-1)
#ani.save("voyager_I_GA_optimisation.gif", writer="pillow", fps=10, dpi=200) # Save as GIF
<Figure size 640x480 with 0 Axes>
In [169]:
# Convert animation to HTML
html_video = ani.to_jshtml()
# Display animation inline
display(HTML(html_video))
In [170]:
fig, ax = plt.subplots(figsize=(9, 7))
plt.close(fig) # This prevents the initial empty figure from showing
plotter = StaticOrbitPlotter(ax)
# Store previously plotted trajectories
plotted_trajectories = []
# Select a colour map (viridis, plasma, coolwarm, rainbow, etc.)
cmap = plt.cm.rainbow # Change this for different color gradients
# Flag to track if we've plotted the real mission trajectory
real_mission_plotted = False
# Improve grid and ticks
ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis='both', which='major', labelsize=12)
ax.tick_params(axis='both', which='minor', labelsize=10)
plt.tight_layout()
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=len(best_solutions_per_generation)+1,
repeat=False, blit=False)
filename_GA_trajectory_per_gen = output_ani.get_path("voyager_I_GA_trajectory_per_gen.gif")
ani.save(filename_GA_trajectory_per_gen, writer="pillow", fps=10, dpi=200)
<Figure size 640x480 with 0 Axes>
In [171]:
def display_saved_gif(filename):
with open(filename, 'rb') as file:
gif_data = file.read()
encoded = b64encode(gif_data).decode('ascii')
html = f'<img src="data:image/gif;base64,{encoded}">'
return HTML(html)
display(display_saved_gif(filename_GA_trajectory_per_gen))
GA Problem Space per Generation¶
In [173]:
# Function to animate the evolution of the genetic algorithm
def animate_3d_ga_contour(population_per_generation):
fig = plt.figure(figsize=(10, 8))
plt.close(fig) # This prevents the initial empty figure from showing
ax = fig.add_subplot(111, projection='3d')
cmap = plt.cm.viridis_r # Color map for deltaV values
cbar = None # Initialise colorbar variable
def update(frame):
nonlocal cbar
ax.clear()
generation_idx = frame
selected_population = population_per_generation[generation_idx]
# Extract dates and normalise relative to nominal mission dates
departure_dates_GA = np.array([ind[0] for ind in selected_population]) - date_departure_mjd
jupiter_flyby_dates_GA = np.array([ind[1] for ind in selected_population]) - date_flyby_jupiter_mjd
saturn_flyby_dates_GA = np.array([ind[2] for ind in selected_population]) - date_flyby_saturn_mjd
# Compute deltaV for coloring
deltaV_values_GA = np.array([find_deltaV_mission(*ind)[0] for ind in selected_population])
scatter = ax.scatter(departure_dates_GA, jupiter_flyby_dates_GA, saturn_flyby_dates_GA,
c=deltaV_values_GA, cmap=parula_map(True), s=10, alpha=0.3,
norm=plt.Normalize(vmin=min(deltaV_values_GA), vmax=min(deltaV_values_GA) * 1.5))
# Mark the best solution for the current generation
best_individual = selected_population[0] # Best individual of the generation
ax.scatter(best_individual[0] - date_departure_mjd,
best_individual[1] - date_flyby_jupiter_mjd,
best_individual[2] - date_flyby_saturn_mjd,
color='red', s=100, marker='x', label="Best Solution")
ax.scatter([0],
[0],
[0],
color='red', s=100, marker='o', edgecolors='black', linewidth=1.2,
label=f'Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg')
# Set labels with increased padding
ax.set_xlabel('Change in Earth Departure Epoch (days)', labelpad=10)
ax.set_ylabel('Change in Jupiter Flyby Epoch (days)', labelpad=10)
ax.set_zlabel('Change in Saturn Flyby Epoch (days)', labelpad=10)
# Set title with increased padding
ax.set_title(f'3D Visualisation of GA Problem Space (Generation {generation_idx})', pad=20)
# Format the scientific notation to be more readable
for axis in [ax.xaxis, ax.yaxis, ax.zaxis]:
formatter = ScalarFormatter(useOffset=False)
formatter.set_scientific(True)
formatter.set_powerlimits((-3, 3)) # Only use scientific notation for exponents outside this range
axis.set_major_formatter(formatter)
# Position the offset text for better visibility
ax.xaxis.get_offset_text().set_position((0, -0.1))
ax.yaxis.get_offset_text().set_position((0, -0.1))
ax.zaxis.get_offset_text().set_position((0, -0.1))
# Ensure the exponent text is horizontally aligned
ax.xaxis.get_offset_text().set_ha('left')
ax.yaxis.get_offset_text().set_ha('left')
ax.zaxis.get_offset_text().set_ha('left')
ax.legend()
limit = 365
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
# Add colorbar only once
if cbar is None:
cbar = fig.colorbar(scatter, ax=ax, shrink=0.6, pad=0.1)
cbar.set_label('Total $\Delta V$ (km/s)')
ani = animation.FuncAnimation(fig, update, frames=len(population_per_generation), repeat=False)
plt.show()
return ani
# Example Usage
ani = animate_3d_ga_contour(population_per_generation)
In [174]:
# Convert animation to HTML
html_video = ani.to_jshtml()
# Display animation inline
display(HTML(html_video))